# Hierarchical Bayesian Modeling for Time-Dependent Inverse Uncertainty
Quantification

Chen Wang11footnotemark: 1

###### Abstract

This paper introduces a novel hierarchical Bayesian model specifically
designed to address challenges in Inverse Uncertainty Quantification (IUQ) for
time-dependent problems in nuclear Thermal Hydraulics (TH) systems. The unique
characteristics of time-dependent data, such as high dimensionality and
correlation in model outputs requires special attention in the IUQ process. By
integrating Gaussian Processes (GP) with Principal Component Analysis (PCA),
we efficiently construct surrogate models that effectively handle the
complexity of dynamic TH systems. Additionally, we incorporate Neural Network
(NN) models for time series regression, enhancing the computational accuracy
and facilitating derivative calculations for efficient posterior sampling
using the Hamiltonian Monte Carlo Method - No U-Turn Sampler (NUTS).

We demonstrate the effectiveness of this hierarchical Bayesian approach using
the transient experiments in the PSBT benchmark. Our results show improved
estimates of Physical Model Parameters’ posterior distributions and a reduced
tendency for over-fitting, compared to conventional single-level Bayesian
models. This approach offers a promising framework for extending IUQ to more
complex, time-dependent problems.

###### keywords:

Hierarchical Bayesian Model , Time-Dependent Modeling , Inverse Uncertainty
Quantification , Thermal Hydraulics , Markov Chain Monte Carlo , Surrogate
Models

[label1]organization=Department of Nuclear, Plasma and Radiological
Engineering, University of Illinois Urbana Champaign

##  1 Introduction

In recent years, the field of nuclear Thermal Hydraulics (TH) has increasingly
relied on computational simulations to predict the behavior of nuclear reactor
systems under various conditions. These simulations are crucial for ensuring
the safety and efficiency of nuclear reactors. However, they are inherently
uncertain due to limitations in accurately modeling complex physical
phenomena. This uncertainty becomes particularly pronounced in time-dependent
simulations, where the dynamics of the system introduce additional layers of
complexity. Inverse Uncertainty Quantification (IUQ) has emerged as a pivotal
methodology to address these uncertainties, specifically by quantifying the
uncertainties associated with Physical Model Parameters (PMPs) in these
systems.

The final report of the BEMUSE project (Perez et al. (2011)) underscores the
necessity of input uncertainty quantification as a pivotal element in
probabilistic uncertainty analysis within Best Estimate plus Uncertainty
(BEPU) frameworks. The prevalent approach in BEPU methodologies involves the
transmission of uncertainties from inputs to model predictions, making the
determination and validation of uncertainty ranges for each parameter vital.
Traditionally, these ranges have been subjectively set based on expert
opinions, highlighting the need for a more empirical and structured method.
This issue was a key focus of the international PREMIUM project (Skorek et al.
(2019)), though it concluded without a unified agreement on the procedures and
guidelines for input uncertainty quantification (UQ). The absence of a
consensus led to the initiation of the OECD/NEA SAPIUM project (Systematic
APproach for Input Uncertainty quantification Methodology) (Baccou et al.
(2019) Baccou et al. (2020)), which sought to establish a methodical approach
to input UQ in nuclear thermal-hydraulics (TH) codes. Following this, in 2022,
the OECD/NEA launched the ATRIUM project (Application Tests for Realization of
Inverse Uncertainty quantification and validation Methodologies in thermal-
hydraulics) Ghione et al. (2023), with the main objective being the practical
implementation of IUQ exercises to demonstrate the efficacy of the SAPIUM
approach.

The landscape of IUQ has been considerably shaped by the progressive
advancements in machine learning (ML) and artificial intelligence (AI). These
technologies have significantly improved the precision and dependability of
simulation models in various sectors, tackling a multitude of intricate
challenges Wang et al. (2023c). Their successful deployment spans several
disciplines, including healthcare (Dong et al. (2021a, b); Chen et al.
(2019)), agriculture (Wu et al. (2022)), transportation (Ma et al. (2022);
Meng et al. (2022); Li et al. (2023a)), clinical research (Xue et al. (2021,
2022)), signal processing (Hu et al. (2022); Li et al. (2023b); Liu et al.
(2021)), structural health monitoring (Liu and Bao (2023)), and industrial
engineering (Li et al. (2023c); Chen et al. (2018, 2023)). The adaptability
and success of ML/AI approaches in these areas provide compelling evidence of
their potential and impart crucial learnings for advancing IUQ initiatives
within the nuclear power industry.

Traditional IUQ approaches, predominantly using single-level Bayesian models,
have provided significant insights into steady-state TH systems. However, they
exhibit limitations when applied to applications with various experimental
conditions and large datasets. Key challenges include handling the high
variability of PMPs under dynamic experimental conditions and avoiding over-
fitting due to unknown model discrepancies or outliers. Recent advancements in
hierarchical Bayesian models Wang et al. (2023a, b) have shown promise in
addressing these issues in steady-state nuclear TH systems, but their
applicability for time-dependent problems remains underexplored.

This study aims to fill this gap by introducing a novel hierarchical Bayesian
model tailored for time-dependent IUQ in nuclear TH systems. The main
objectives are twofold: firstly, to develop a model capable of handling the
unique challenges of time-dependent data, such as high dimensionality and
correlation; and secondly, to demonstrate the hierarchical model’s efficacy in
providing more accurate and reliable estimates of PMPs’ posterior
distributions, while reducing the tendency for over-fitting. The study
leverages advanced statistical techniques, including Gaussian Processes (GP)
integrated with Principal Component Analysis (PCA), and Neural Network (NN)
models for time series data regression. Hamiltonian Monte Carlo Method - No
U-Turn Sampler (NUTS) are used for efficient posterior sampling. The method’s
performance is evaluated using the transient experiments in the PSBT
benchmark, setting a new precedent in the application of hierarchical Bayesian
methods to time-dependent IUQ in nuclear TH systems.

The rest of the paper is organized as follows. Section 2 will give an overview
of Hierarchical Bayesian model for IUQ. Section 3 will introduce the overview
of the PSBT benchmark and TRACE modeling. Section 4 will introduce the
surrogate models for time-dependent problems, and then in Section 5, the
hierarchical framework is applied to a case study for TRACE physical model
parameters using the PSBT benchmark data. Section 6 will be the summary.

##  2 Hierarchical Bayesian Framework for Inverse Uncertainty Quantification

Hierarchical model is also called a multi-level model or mixed-effect model. A
simple illustration of the structure of the hierarchical model is shown in
Figure 1. Observations $y$ are from different clusters determined by (1)
cluster parameter $c_{m}$ which might be different among clusters, and (2)
parameter $\bm{\theta}$ shared across clusters. If the probability
distribution of $c_{m}$ can be parameterized by $\bm{\Sigma}_{c}$,
$\bm{\Sigma}_{c}$ can be seen as the cluster-specific parameter, or per-group
parameter, to be estimated (Wang (2020)).

Figure 1: Non-hierarchical and hierarchical models

Consider a simple example for Figure 1, where the observations are normally
distributed. In this case, the parameters we need to estimate are the mean and
the variance. Suppose they have the same variance $\sigma^{2}$ for all
clusters but different means $\mu_{i}$ for each cluster. So the $\sigma^{2}$
is a shared parameter ($\theta$ in Figure 1) and $\mu_{i}$ is per-group
parameter ($c_{i}$ in Figure 1). Now we would like to specify the distribution
over the cluster-specific $\mu_{i}$. We can assume the distribution also
normal, or other distributions depending on the actual scenario. If normal
distribution is assumed, two parameters (global mean $\mu$ and global standard
deviation $\sigma_{y}$, correspond to the $\Sigma_{c}$ in Figure 1) would be
required. So we can describe parameter $\mu_{i}$ and observations $y_{ij}$ as:

$\mu_{i}\sim N(\mu,\sigma_{y}^{2})$ $y_{ij}\sim N(\mu_{i},\sigma^{2})$

Now we can re-arrange the structure of the above hierarchical model to better
illustrate the process of the above two equations: the per-group parameters
are generated from shared parameters, and observations are generated from per-
group parameters, as shown in Figure 2.

Figure 2: Illustration of a classical hierarchical model structure

In summary, the hierarchical model comes in handy when we know that some
cluster-level attributes exist and the structure of the data can be reasonably
modeled or assumed. When data follows a hierarchical structure, a simple non-
hierarchical model would be inappropriate because (1) a model with few
parameters cannot fit the dataset accurately, and (2) a model with many
parameters will “overfit” the dataset. In contrast, a hierarchical model can
have enough parameters to fit the data, while has a population distribution to
model the dependence of the parameters. Thus “overfitting” can be avoided
because the parameters are “constrained“ by the population distribution so
they cannot fit the data exactly.

The hierarchical model structure in Figure 1 can be plotted into a more
succinct and formal graphical representation, as shown in Figure 3. In the
figure, $p(\theta)$ is the prior distribution of $\theta$, $i$ is the cluster
index, $y$ is the observation, and there are $N$ observations in total. the
per-group parameter $b_{i}$ is governed by its distribution $\Sigma_{b}$,
which has prior distribution $p(\Sigma_{b})$.

Figure 3: Graphical representation of a hierarchical model

It should be noted that in this model, our primary goal is to learn about
$\theta$ and $\Sigma_{b}$, not the individual cluster-specific parameters
$b_{i}$, but $b_{i}$ also has to be estimated in order to estimate
$\Sigma_{b}$. So we need to marginalize over the $b_{i}$ parameters in after
we have obtained the joint posterior distribution. There are generally two
approaches to estimate the posterior, we can either use MLE for point
estimation or use Bayesian inference to get full posterior distributions for
these parameters.

In the MLE method, the likelihood of $\theta$ and $\Sigma_{b}$, marginal over
$b$, can be expressed as:

$L(\bm{\Sigma}_{b},\bm{\theta};\bm{y})=\int_{b}P(\bm{y}\mid\bm{\theta},\bm{b},i)P(\bm{b}|\bm{\Sigma}_{b})\,db$
(1)

If full posterior distribution is desired, we can incorporate the prior
distribution and use Bayes’ rule. The prior distribution can be decomposed
into the following form:

$P(\bm{\Sigma}_{b},\bm{\theta},b_{i})=P(b_{i}\mid\bm{\Sigma}_{b})\cdot
P(\bm{\Sigma}_{b},\bm{\theta})$ (2)

which is based on the fact that $\Sigma_{b}$ influences observations only
through $b_{i}$. So the marginalized joint posterior distribution can be
expressed as:

$P(\bm{\Sigma}_{b},\bm{\theta}\mid\bm{y})\propto\int_{b}P(\bm{y}\mid\bm{\theta},\bm{b},i)P(\bm{b}\mid\bm{\Sigma}_{b})\cdot
P(\bm{\theta})P(\bm{\Sigma}_{b})\,db$ (3)

Hierarchical models provide a flexible framework for modeling the complex
interactions but come with higher computational requirements in inference. The
parameter space to be estimated can be high-dimensional: we need to estimate
cluster-specific parameter $b_{i}$, which may contain up to hundreds.
Traditional random-walk based MCMC algorithms such as MH or advanced MH
algorithms tailored for relatively high dimensional problems still suffer from
serious convergence issues because that the random behavior of the proposal
function in very inefficient in high-dimension domains. HMC and NUTS would be
useful in such high-dimensional conditions.

##  3 Transient Experiments in PSBT Benchmark

OECD/NRC benchmark based on NUPEC PWR subchannel and bundle tests (PSBT) is
designed for validation purposesof void distribution in subchannel and PWR
bundle and pre-diction of departure form nucleate boiling Uncertainty
Quantification for Steady-State PSBT Benchmark using Surrogate Models. The
void distribution benchmark in PSBT includes transient bundle benchmark, which
can be applied to system TH codes to assess their capabilities of predicting
the void generation during transients (Rubin et al. (2012)). The experimental
data in these transients include X-ray densitometer measurements of void
fraction (chordal averaged) at three axial elevations. The averaging is over
the four central subchannels. Data is collected for four transient scenarios:
Power Increase (PI), Flow Reduction (FR), depressurization (DP), Temperature
Increase (TI), and at three different assembly types 5, 6, and 7. All 5, 6,
and 7 assemblies are $5\times 5$ rob bundles while 5 and 6 have typical cells
and 7 has thimble cells. They also have different axial and radial power
distributions. The PSBT benchmark has been used extensively in many IUQ
related applications (Borowiec et al. (2017); Wang et al. (2018a))

In this work, we managed to build TRACE models for power increase, flow
reduction, and temperature increase transients for all three assemblies.
Figure 4 shows the variation of boundary conditions in 7T-flow reduction
transient: during the flow reduction transient, flow rate was reduced and
other boundary conditions were maintained at the same level.

Figure 4: Boundary conditions during the 7T-flow reduction transient

Figure 5, 6, and 7 show comparisons between TRACE predicted void fraction and
measured void fraction for the three assemblies, respectively. Dots represent
experiment measurements and lines are TRACE simulations. We can see that in
some cases, TRACE performs well and have good agreement with experimental
data, while there are also obvious model discrepancies in some transients.

Figure 5: Transients simulation in Assembly 5. (Dots: Experiment, Lines:
Simulation) Figure 6: Transients simulation in Assembly 6. (Dots: Experiment,
Lines: Simulation) Figure 7: Transients simulation in Assembly 7. (Dots:
Experiment, Lines: Simulation)

Since the transient and the steady-state experiments share same assemblies and
similar boundary conditions, it is reasonable to use the same physical model
parameters as we have selected in previous similar studies (Wang et al.
(2017a, 2019a, 2018b, 2019b)). The four parameters are listed in Table 1.

Table 1: List of 4 selected physical model parameters in TRACE Parameter Number | Definition
---|---
P1008 | Single phase liquid to wall heat transfer coefficient
P1012 | Subcooled boiling heat transfer coefficient
P1022 | Wall drag coefficient
P1028 | Interfacial drag (bubbly/slug Rod Bundle-Bestion) coefficient

##  4 Surrogate Model for Time-Dependent Thermal-Hydraulics Systems

Surrogate model, metamodel or emulator plays a vital role in the Bayesian
calibration of computer models when the computer models are computationally
expensive, such as nuclear reactor simulations. The surrogate models are built
from limited number of runs of the full model at specifically selected values
of the input parameters and a learning algorithm. The Design of Experiment
(DoE) or computer experiments design Santner et al. (2013), in the context of
surrogate models, aims to select a structured set of tests of the computer
model, such that a good tradeoff between the accuracy of the surrogate model
and the number of tests can be achieved. A good coverage of the input space
will ensure the accuracy of the emulator while too many samples will cause
extra computational burden, which is against the purpose of using the
surrogate models. Latin Hypercube Sampling (LHS) is a widely-used method to
generate samples from a multidimensional distribution. It has good space-
coverage property and will be used in this work.

Many types of surrogate models have been utilized in the Bayesian calibration
process, such as the Polynomial Chaos Expansion (PCE) and Stochastic
Collocation (SC) used by Wu et al. (2017a) Wu and Kozlowski (2017), the
Multivariate Adaptive Regression Splines (MARS) proposed by Stripling et al.
(2013), the Function Factorization with Gaussian Process (FFGP) developed by
Yurko et al. (2015), and the deep Neural Networks by Liu et al. (2018). There
are many more applications that cannot be included here. Different types of
surrogate models have different characteristics, so they may fit in different
engineering scenarios. Some factors to consider in choosing a proper one
include the dimension of inputs, non-linearity of the input-output
relationship, and time-dependent or steady-state problem, etc.

###  4.1 Gaussian Processes

Gaussian Processes have been widely used as surrogate models since the seminal
work of Bayesian calibration by Kennedy and O’Hagan (2001) utilized GP
modeling. In this section, we will firstly go through the theoretical
backgrounds of GP ( Wang et al. (2019a, 2017b)), and look into several
characteristics of GP.

A Gaussian process can be used to describe a distribution over functions. It
is normally defined as a collection of random variables, any finite number of
which have a joint Gaussian distribution ( Rasmussen (2004)). A GP is
completely defined by its mean function $m(\bm{x})$ and covariance function
$k(\bm{x},\bm{x}^{\prime})$, and can be written as:

$f(\bm{x})\sim\mathcal{GP}(m(\bm{x}),k(\bm{x},\bm{x}^{\prime}))$ (4)

where the random variables in the process represent the value of the function
$f(\bm{x})$ at location $\bm{x}$. For notation simplicity we will treat the
mean function as zero at this time. Later we will see this is not a drastic
limitation because the posterior mean is not confined to be zero, and the
prior mean function does not have major influence in interpolation.

Given training inputs and outputs $\left\\{(\bm{x_{i}},f_{i})\mid
i=1,...n\right\\}$, we want to predict the test outputs $\mathbf{f_{*}}$ in
test points $X_{*}$. The prediction outputs (including training outputs
$\mathbf{f}$ and test outputs $\mathbf{f_{*}}$) are jointly distributed as a
multivariate normal distribution:

$\begin{bmatrix}\mathbf{f}\\\ \mathbf{f}_{*}\end{bmatrix}\sim
N\Bigg{(}0,\begin{bmatrix}K(X,X),K(X,X_{*})\\\
K(X_{*},X),K(X_{*},X_{*})\end{bmatrix}\Bigg{)}$ (5)

where $K(X,X_{*})$ denotes the matrix of the covariances evaluated at all
pairs of training and test points, and similar definitions apply to $K(X,X)$,
$K(X_{*},X_{*})$, and $K(X_{*},X)$. We are interested in the conditional
probability of $\\{\mathbf{f}_{*}|X_{*},X,\mathbf{f}\\}$, which, fortunately,
is simple to obtain due to properties of multivariate normal distributions.

$\begin{split}\mathbf{f}_{*}\mid X_{*},X,\mathbf{f}\sim
N\Big{(}&K(X_{*},X)K(X,X)^{-1}\mathbf{f},\\\
&K(X_{*},X_{*})-K(X_{*},X)K(X,X)^{-1}K(X,X_{*})\Big{)}\end{split}$ (6)

Now we can evaluate the function values $\mathbf{f}_{*}$ according to Equation
(6) by calculating the corresponding mean and the covariance matrix. The
posterior means are used as prediction values. The covariance matrix represent
the level of confidence in the prediction, and it will enter the covariance
term $\bm{\Sigma}$ in Equation LABEL:eqa:sigma_total, so that the covariance
term in the likelihood will include the error information from both the
experiment and the emulator:
$\bm{\Sigma}=\bm{\Sigma_{\text{exp}}}+\bm{\Sigma_{\text{emu}}}$. This
treatment is important because inaccurate predictions by the emulator can be
seen as assigned a much smaller weight in evaluating the likelihood in
Equation LABEL:eqa:sigma_total. This unique characteristic of GP makes it a
preferable choice of an surrogate model in the Bayesian Calibration process.

The covariance function is a crucial ingredient in a GP predictor as it
encodes the information of nearness or similarity between data points. This
information is crucial because in regression tasks we normally assume closer
input points are more likely to have similar target values, and thus the
training points which are nearer to a test point should have more influence on
the prediction at that point. The choice of the covariance function should
depend on the assumption of the function which we wish to learn. In most of
the cases and in this work it is defined as a function of the Euclidean
distance between two points, with several hyperparameters. For example, power-
exponential covariance function:

$\text{cov}(f(\bm{x}_{p}),f(\bm{x}_{q}))=k(\bm{x}_{p},\bm{x}_{q})=\sigma^{2}\cdot\exp(-\frac{1}{2l^{2}}\mid\bm{x}_{p}-\bm{x}_{q}\mid^{r})$
(7)

The hyperparameters $\sigma,l,$ and $r$ define the proprieties of the
covariance function such as magnitude, shape, and smoothness. They should be
properly estimated to obtain the best possible prediction performance.

There is also Gaussian Kernel, where $r$ is simply changed to 2:

$k(\bm{x}_{p},\bm{x}_{q})=\sigma^{2}\cdot\exp(-\frac{1}{2l^{2}}\mid\bm{x}_{p}-\bm{x}_{q}\mid^{2})$
(8)

The Matérn class of kernel is given by:

$k(\bm{x}_{p},\bm{x}_{q})=\frac{1-\nu}{\Gamma(\nu)}(\frac{\sqrt{2\nu}\mid\bm{x}_{p}-\bm{x}_{q}\mid}{l})^{\nu}K_{\nu}(\frac{\sqrt{2\nu}\mid\bm{x}_{p}-\bm{x}_{q}\mid}{l})$
(9)

with positive parameters $\nu$ and $l$, and $K_{\nu}$ is a modified Bessel
function (Rasmussen (2004)). $\nu=3/2$ and $\nu=5/2$ are popular choices.
Different kernels may fit in different scenarios.

Kernels functions play important roles in GP modeling. A GP model with the
same parameter estimation method but different kernel functions may have
different results. This can be seen in Figure 8.

Figure 8: GP regression with RBF and Matérn ($\nu=5/2$) kernels

A simple function $f(x)=sin(x)*x$ is created and several points on $x$ are
sampled as training data (red circles). As mentioned previously, the GP model
yields both the mean and the standard deviation of the response, which allows
us to plot the $95\%$ CI for each unobserved data point, see blue shaded area
in Figure 8. Despite the prediction results (blue line) by two GP models
showed good performance compared the the observations, the standard deviations
yielded by two models are quite different. The RBF kernel showed small
standard deviation while the model with Matérn kernel has relatively large
standard deviation. So extra caution should be paid when incorporating the
error information of GP model to the total covariance $\bm{\Sigma}$ in
Bayesian calibration. It is reasonably that different models have different
confidence in their prediction, but the difference in the resulting total
covariance matrix may lead to very different posterior distributions. So it is
recommended that the uncertainty brought by GP surrogate should not be
explicitly included in the total covariance.

The GP model in this work is coded by the scikit-learn (Buitinck et al.
(2013)) open-source package in Python. The default parameter estimation method
in scikit-learn is “L-BGFS-B” from scipy (Virtanen et al. (2020)). The
parameter estimation aims at maximizing the log-marginal-likelihood (LML). The
LML of GP may have multiple local optima, so the algorithm will restart
automatically for a certain number of times to ensure an optimized set of
hyper-parameter can be found. Cross Validation (CV) can also be used to
estimate the hyper-parameters. Details about various parameter estimation in
GP can be found in (Rasmussen (2004)). The accuracy of the GP emulator needs
to be assessed before use. We are interested in the predictive accuracy at
untried points, which can be done by quantifying the predictive error at an
additional set of validation data (Wu et al. (2017c, b)).

###  4.2 Principal Component Analysis and Gaussian Processes

When dealing with the time series regression problems, traditional regression
methods becomes challenging as the highly-correlated output usually needs
special treatment by many separate GP models which increases computational
burden significantly. A multivariate GP surrogate model is required in this
case instead of independent model for each output variable. For example, if
there are a hundred measurements in a series, then all the hundred outputs
need to be processed by the surrogate model, requiring significant higher
training time. A dimension reduction of the high-dimensional output helps
reduce the total computational burden and maintain the correlations of
outputs. In this section, foundations of the PCA will be firstly introduced. A
combined time series regression model by GP and PCA will be explained.

PCA can be viewed as a projection method which projects high-dimension data
into a lower-dimension subspace such that the projected data are uncorrelated
and have maximized explained variance ( Wu et al. (2018)Wang et al. (2018b)).
PCA can be done through Eigen-decomposition of the covariance matrix or
Singular Value Decomposition (SVD) of the data matrix.

A brief mathematical introduction to PCA will be given here. Suppose our data
matrix $\bm{A}$ has the shape $p\times N$, where $N$ stands for number of
samples, $p$ stands for number of features (dimension of each output). This
data matrix can be seen as the training dataset for the surrogate model.

PCA is a linear transformation method where we transform the data matrix
$\bm{A}$ by

$\bm{P}\bm{A}=\bm{B}$ (10)

where $\bm{P}$ is a transformation matrix and $\bm{B}$ is a transformed data
matrix. The goal of PCA is to find the proper $\bm{P}$ that makes rows of
$\bm{B}$ uncorrelated and this can be done by SVD. If A is already centered
(subtract row mean for each element), the SVD of A can be given by:

$\bm{A}=\bm{U}\bm{\Lambda}\bm{V}^{T}$ (11)

where $\bm{U}$ is an orthogonal matrix whose columns are the left-singular
vectors of $\bm{A}$, $\bm{V}$ is an orthogonal matrix whose columns are the
right-singular vectors of $bmA$, and $\bm{\Lambda}$ is a $p\times N$ diagonal
matrix whose diagonal entries are the singular values of $\bm{A}$. The
magnitudes of diagonal entries of $\bm{\Lambda}$ are then arranged in
descending order and large singular values corresponds to important features
in $\bm{A}$ (Shlens (2014)). Choose $\bm{P}$ as $\bm{U}^{T}$, so that

$\bm{P}\bm{A}=\bm{\Lambda}\bm{V}^{T}=\bm{B}$ (12)

The covariance matrix of $\bm{B}$ is diagonal. The rows of $\bm{P}$ and the
rows of $\bm{B}$ are called Principal Components (PC) and PC scores,
respectively.

The dimension reduction can then be performed by selecting a proper dimension
of the principal subspace $p^{*}$, which is much smaller than $p$, but is able
to make the principal components explain most of the variances. The first
$p^{*}$ PCs form a $p^{*}\times p$ transformation matrix $\bm{P}^{*}$:

$\bm{P}^{*}\bm{A}=\bm{B}^{*}$ (13)

where $\bm{B}^{*}$ is a $p^{*}\times N$ matrix whose rows represent new
variables after dimension reduction. So, matrix $\bm{B}^{*}$ has much lower
dimension than the original data matrix, and can be used as the training data
for GP.

The left matrix multiplication of the data matrix can be seen as rotating,
moving, and scaling the coordinate system of the original data. Suppose we are
looking at a 2 dimensional input shown in Figure 9, and we want to conduct
dimension reduction on the data. It can be seen that the data has the
potential to be reduced to one dimension because $x$ and $y$ are linearly
related with some noise, knowing one of them can let us know roughly about the
other. So the axes can be moved, rotated, and scaled as the red arrows show in
Figure 9. The longest arrow (axis) can explain most of the variance of the
data, while the data variation along the second arrow (axis) is small compared
to the first.

Figure 9: Illustration of PCA in 2D space Figure 10: Comparison of original
data and transformed data by PCA using one PC

The transformed data $\bm{B}$ obtained by PCA can be seen in the lower plot of
the Figure 9. The PCs are actually the projection of each data point onto the
principal axes. From the upper figure we can see that PC1 explains most of the
variation so we select PC1 only and transform it back to the original input
space, the results are shown in Figure 10. We can see the linear relationship
between $x$ and $y$ is still captured after dropping PC2, and the transformed
data using only 1 principal component is “good enough” to represent the
original data.

From Figure 10, we can have a better understanding of that PCA does for
dimension reduction: the information along the least important axes is
removed, and the component(s) of data which explains most of the variance
is(are) left. The fraction of the variance that is cut out can be seen as a
measure of how much information is discarded by the dimension reduction of
PCA. This metric is usually used to determine the number of PCs to select.

In the reduced dimension, we would be able to combine GP model with PCA to
construct a regression model. The flowchart of using GP and PCA as regression
model for the high-dimension and high-correlation problem is shown in Figure
11.

Figure 11: Flowchart of the GP + PCA regression model for high-dimension and
high-correlation outputs

The dataset should be split into a training set and validation set so that
over-fitting can be avoided and the model performance can be objectively
quantified. In the training set with $N$ samples, suppose the the input and
the output has the dimension d_in and d_out, respectively. The output we are
dealing here is time-dependent data so it has high dimensions and high
correlation. Applying PCA to the output and selecting proper number of
principal components $p^{*}$, we can reduce the output dimension to $p^{*}$.
The selection of number of PCs should be based on the fraction of variance
explained by the PCs. A GP model can then be built here based on the original
input (d_in, N) and the output ($p^{*}$, N) in the reduced dimension.

In regression, we are interested in predicting the output at an new untried
point, and the prediction by GP model here gives us the output in the reduced
dimension with shape ($p^{*},1$). Inverse of PCA can then be applied here to
transform the data into the original output space. The inverse can be done
simply by left multiplication of $\bm{P}^{-1}$ to Equation 12. If the
predicted output is compared with the new output in the validation set, the
performance of the overall model can then be quantified.

In summary, the combined GP and PCA model requires building one PCA model and
$p^{*}$ GP models, which reduces the computational burden of building d_out GP
models. The dimension reduction becomes possible because these time-dependent
outputs are highly correlated.

###  4.3 Artificial Neural Network

Artificial Neural Network (ANN) is a powerful tool for regression tasks for
its capability to model complex patterns in datasets using multiple hidden
layers and non-linear activation functions. In many cases, Neural Network may
seem like an overkill as regression models since it is mostly known for its
capability of handling very complex and computationally expensive jobs such as
image processing, natural language process, artificial intelligence, etc.
However, NN or even Deep NN is nothing but repeated combinations of simple
linear and non-linear transformations. Their structures can be designed in a
very sophisticated way but the smallest unit in DNN is still a perceptron. So
ANN models can be seen as convenient and powerful regression models if
properly tuned.

Perceptron (or neuron) is the smallest unit in a neural network. The output of
a perceptron is simply an activation function applied to linear combination of
inputs. Say if the input $\bm{x}$ has $n$ dimensions, now the output from a
perceptron if $f(b+x_{1}w_{1}+x_{2}w_{2}+...+x_{n}w_{n})$, where $b,w_{i}$ are
weights and $f(.)$ is the activation function. Many activation functions are
available and they have different characteristics to be used in different
scenarios. Most commonly used activation function for regression purposes
include linear or Rectified Linear Units (ReLU) (Ramachandran et al. (2017)).
Linear activation provides simple unbounded output and is easy to solve, but
is limited in its complexity and have less power in complex problems. ReLU is
widely used in many ANN structures and it ranges from 0 to infinity. It has
the form:

$f(x)=\begin{cases}0,x<0\\\ x,x\geq 0\end{cases}$ (14)

which gives an output itself if positive and 0 otherwise. It is a good
approximator and any function can be approximated with a combination of Relu
Ramachandran et al. (2017).

Mathematically, a one layer network has the form:

$\bm{y}=\bm{w}^{T}\bm{x}$ (15)

and a two layer network has the form:

$\bm{y}=\bm{w}^{T}\sigma_{1}(\bm{A}_{1}\bm{x}+\bm{b}_{1})$ (16)

where $\bm{A}$ is the weight matrix and $\bm{b}$ is the intercept vector of
that layer. Iterating the process, a fully connected network with $L$ layers
has the form:

$\bm{y}=\bm{w}^{T}\sigma_{1}(\bm{A}_{1}\sigma_{2}(...\bm{A}_{L-2}\sigma_{L-1}(\bm{A}_{L-1}\bm{x}+\bm{b}_{L-1})+\bm{b}_{L-2}...)+\bm{b}_{1})$
(17)

The objective function of NN model in this regression task can be defined as
the sum of squared errors: $\Sigma(\hat{y}_{i}-y_{true})^{2}$, which is
usually referred to as L2 norm loss function. The parameters in the model
($\bm{A},\bm{b}$) can then be estimated by minimizing the objective function
with respective to $\bm{A}$ and $\bm{b}$. Taking derivatives of equation
LABEL:fig:ann_1 with respect to these weight parameters is straightforward and
can be scaled to any finite depth because of the chain rule. The method of
adjusting weight parameters to minimize the objective function using the chain
rule is called backpropagation (Hecht-Nielsen (1992)). With the development of
computer science and computational power, the backpropagation algorithm can be
applied to very large-scale neural networks and has achieved successes in many
fields.

ANN model should be used with special caution in regression tasks because it
is very prone to overfitting due to its complex structure. To avoid
overfitting, there are generally two methods to considered when the number of
available data points can not be increased:

  * 1.

Weight Regularization. Add a penalty term to the error function during the
training process, so that weight parameters can be constrained. Because large
weighs tend to cause sharp transitions in the activation functions and thus
large changes in output for small changes in inputs (Reed and MarksII (1999)).

  * 2.

Add dropout layers. During the training process, we randomly drop units (along
with their connections) from the neural network so that it prevents units from
adapting too much with each other (Srivastava et al. (2014)).

Both two methods have been widely used to prevent overfitting. One additional
benefit of adding dropout layers is that model uncertainty can also be
calculated besides the point estimate (Gal and Ghahramani (2016)), which gives
us confidence interval for the prediction. So dropout method is preferred in
many application where uncertainty information is needed such as reinforcement
learning and Bayesian optimization (Snoek et al. (2012)).

In this work, there are two main motivations of using ANN as surrogate models.
The first is that we need to deal with high dimensional output in the time-
series data, and ANN is capable of dealing with high dimensions with its
flexible structure. The second is that some efficient MCMC algorithms require
derivatives of the model prediction, which is easy to obtain from Equation 17.
The ANN model is coded using the TensorFlow package (Abadi et al. (2015)) in
Python.

###  4.4 Construction and Validation of Surrogate Model

####  4.4.1 GP and PCA Model

A surrogate model needs to be built for each transient. In this section, we
take 6T-flow reduction transient as an example to illustrate how the surrogate
model is built and validated. In previous section, we have explained in detail
the process of building the GP + PCA model for the time series regression
model. Now we will go through the process using the 6T-FR transient as an
example.

In the 6T-FR transient, there are 40 time steps and 3 outputs, so the output
dimension is $3\times 40=120$. Now we firstly use PCA as a dimension reduction
tool to reduce the original dimension (120) to a much smaller number, and this
number is just the number of principal components we will select. 400 samples
are firstly drawn using LHS in the range of $\theta\in(0,5)$ as the training
set. PCA is then conducted on this $400\times 120$ matrix and the ratio of
explained variance versus the number of principal components is shown in
Figure 12. We can see that 4 principle components can lead to a very good
0.999 percentage of the variance, which means that the selected 4 PCs are a
good enough representation of the original 120 dimensional data.

Figure 12: Ratio of explained variance v.s. number of principal components

To further analyze the effect of the number of PCs, we can transform the
reduced data back to the original space, and see if there is any discrepancy
between the transformed data and the original data. Mathematically, this is
equivalent to comparing $\bm{P}^{-1}\bm{B}^{*}$ with original data matrix
$\bm{A}$. (Recall that in PCA, $\bm{A}=\bm{P}^{-1}B$, $\bm{P}^{*}A=\bm{B}^{*}$
and $\bm{B}^{*}$ is the matrix after dimension reduction. See equation 10.)
The comparison results for one of the samples are shown in Figure 13, where 2,
3, and 4 PCs are tested. In the upper plot of Figure 13, we can see the red
squares are slightly different than the dotted lines, indicating a loss of
accuracy using 2 PCs. The results with 3PCs show better agreement but slight
difference can still be observed. When 4 PCs are used, we can see that the
transformed data and the original date agree almost perfectly.

The same analysis needs to be conducted for the other 8 transients because
they have different profiles and may need more PCs to achieve good accuracy.
Results will not be repeated here for simplicity. Finally, we found out that 5
PCs are enough for very good accuracy for all transients, and we will use 5
PCs as the output to train GP regression models. We have successfully reduced
the output dimension from 120 to 5 using PCA thanks to the correlations in the
time-dependent data.

Figure 13: Comparison between original data (TRACE) and back-transformed data
using 2, 3, and 4 PCs. The discrepancy between the curves are decreasing as
more PCs are used. Figure 14: Convergence study of GP PCA model

Using the 5 PCs, we can train GP models whose inputs are physical models
parameters and outputs are these PCs. Now we can quantify the accuracy of the
combined GP + PCA model as a whole. 50 randomly selected samples are used as a
testing set. The error is calculated as the difference between original TRACE
simulation results and the predicted results by GP + PCA model, on the testing
set. Note that each sample corresponds to a time series so the error is
calculated for multiple time points in a series. The mean absolute error
versus the number of training samples is plotted in Figure 14. We can see that
after 200 samples, the errors are below 1% and the curve remains relatively
flat. But 400 samples lead to the lowest MAE so we will use 400 samples for
this study. This may seem computationally expensive to draw so many samples
from TRACE, however, since this is a one-time expense and the transients
simulations don’t take a very long time so it is acceptable to sacrifice
several hours to achieve better accuracy.

###  4.5 Artificial Neural Network Model

Although GP+PCA model showed very good predictive performance on the high-
dimension high-correlation output, it has two disadvantages: (1) the model
structure is complex and not friendly to users; (2) it is not easy to
calculate the derivative of the model responses with respective to inputs,
which is required in NUTS algorithm for efficient posterior sampling. In this
work, we developed an ANN model, which is capable of handling high-dimension
output, has a straightforward model structure, and provides convenient
derivative information.

The downside of ANN is that its model structure is so flexible that it
requires a lot of tuning work in order to find an optimal structure for the
given task. The current problem requires a 4-dimensional input and
120-dimensional output (for 7T-FR transient only, output shape may be
different for other transients). So we employ a fully connected NN with 4
nodes in the first(input) layer and 120 nodes in the last(output) layer and
use linear activation function in the output layer since it is a regression
problem. The structure of the hidden layers are determined by several trials
and errors, the structure that leads to the best overall prediction
performance in the validation set is selected. The final model has two hidden
layers with 32 nodes in each, so the final NN model is connected by $4\times
32\times 32\times 120$ nodes. L2 norm is applied to reduce the overfitting
issue of the NN model.

The model is trained on the same 400 samples used for GP+PCA in the previous
section, extra 50 random samples are used as the validation set. Figure 15
shows the training history of the ANN model, and we can see that both
regression MAE and MSE decreased to a low level after a certain number of
iterations. The error in the validation set (dotted lines) also converged and
overfitting didn’t occur.

Figure 15: Training history of ANN

The performance of the ANN model and GP+PCA model are compared in the
convergence study shown in Figure 16. MAEs of void fraction predictions are
compared using different number of samples for training. We can see that the
GP+PCA model outperforms the ANN model when there are less than 250 samples,
while the ANN model showed better predictive capability with more data.

Figure 16: Convergence study for GP+PCA and ANN

##  5 Results of Hierarchical Bayesian Model for Time-Dependent Problems

###  5.1 Covariance in Time-Dependent Data

Now that we have solved the issue regarding surrogate models for time-
dependent data, let’s go back to the model updating equation to review some
assumptions. In the model updating equation when model discrepancy is not
considered, the experimental observations actually follow multivariate normal
distribution
$\bm{X}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma}),\bm{X}\in\mathcal{R}^{k}$, where
$\bm{\mu}$ can be considered as the true noise-free value provided by
simulation model, and $\bm{\Sigma}$ is the covariance matrix, which contains
uncertainty information of experiment and code. Then the probability density
function of $\bm{X}$ follows:

$f(\bm{x})=(2\pi)^{-\frac{k}{2}}\det(\bm{\Sigma})^{-\frac{1}{2}}\exp(\frac{1}{2}(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu}))$

Traditionally, we have been treating $\bm{\Sigma}$ as a diagonal matrix,
meaning that we assumed $\bm{X}$ is independently normally distributed and
ignored $cov(X_{i},X_{j})$ when $i\neq j$. This is a valid assumption in most
steady-state experiments, since measurements were usually conducted
independently. However, in transient cases, there can be high correlation
between the measurement errors at different time steps, causing $\bm{\Sigma}$
to be non-diagonal. The correlated measurement error has been studied in many
fields such as hydrologic models by Tiedeman and Green (2013)), Meteorology
models by Stewart et al. (2008), etc. There can be various reasons for the
error correlation, such as spatially and temporally correlated anomalies in
the underlying physical phenomena, multiple observations derived from a single
measurement, etc (Tiedeman and Green (2013)). The correlated observation error
in the transient experiments in the PSBT benchmark can be observed from Figure
5, 6, and 7. For example, in the flow reduction transient in Figure 7, the
boundary conditions are stable after around 163$s$, so the void fractions
measured at three locations are expected to be also stable (normally
distributed at certain levels). However, we can see clear serial correlation
in those time series, indicating correlated errors.

The off-diagonal element in covariance matrix is calculated by
$\Sigma_{i,j}=E[(X_{i}-\mu_{i})(X_{j}-\mu_{j})]$, which requires a bunch of
samples of the measured time series data. The samples cannot be obtained
realistically since experiments cannot be repeated, but we can “simulate” the
experiment process and the error by our simulation model TRACE, using the
uncertainty information from boundary conditions. Since the measured
uncertainty in boundary conditions is known to us, we can propagate its
uncertainty through the code and treat the resulting model responses’
uncertainty as a proxy of the real measurement uncertainty in time series.
Since the code is deterministic we also need to add independently distributed
error to it during the calibration process. It should be noted that this
method can be seen as an over-estimation of the covariance matrix, because the
correlation would exist in the whole time series using this method. However,
from what we can observe in Figure 7, the correlations only exist in multiple
time steps and not the whole time series. There are also other methods
available for estimating the observation covariance matrix, but none are
without fault (Waller et al. (2016)Stewart et al. (2013)).

Figure 17 shows the uncertainty of void fraction caused by boundary conditions
uncertainties in 5T-Flow Reduction transient. Using this ensemble of time-
series realizations, we can calculate the covariance matrix to be used in the
posterior distribution.

Figure 17: Simulated uncertainty in 5T-Flow Reduction transient

The heatmap of the covariance matrix calculated by samples of the ’Lower’
output (yellow curves) in Figure 17 is shown in Figure 18. We can see that the
correlations between time steps do exist and we should take that into
consideration. Same procedure are applied to the other 8 transients and
covariance matrix

Figure 18: Heatmap of covariance matrix using data from the Lower measurement
location in Figure 17.

###  5.2 Posterior Sampling Results and Validation

Now we can use the obtained covariance matrix in the formulation of posterior,
and sample posterior distributions using the MCMC algorithm. In this section,
we will demonstrate the calibration results using a single transient 7T-Flow
Reduction, and compare the results when covariance is included and not
included.

Because of the usage of the ANN model, the efficient NUTS algorithm can be
used for posterior sampling. We firstly show results when covariance is not
considered. In this case, the procedure is exactly the same as calibrating
steady-state data. In this transient experiment, 151 data points are
available, so we can simply calibrate the physical model parameters against
these data, with the covariance matrix being diagonal. Figure 19 shows
posterior distributions and trace plots of the 4 physical model parameters,
and Figure 20 shows the pair-wise joint distribution among the four
parameters. We can see that the chain is converged and mixes well. Two
parameters ‘P1008’ and ‘P1022’ tend to 0 which can be an indication of model
discrepancy. ‘P1028’ the wall drag coefficient also deviates far away from its
nominal values. Here the range of ‘P1028’ is larger than the original range of
(0,5) because the iterative re-sampling procedure is used to ensure the
posterior can be fully included. Correlation of input parameters exist between
‘P1008’ and ‘P1012’, and between ‘P1012’ and ‘1028’, as can be seen in Figure
20.

Now we include the covariance matrix in the posterior sampling, the posterior
distributions, trace plots and pair-wise joint distributions are shown in
Figure 21 and 22. We can see the results are very different from Figure 19 and
20, where covariance is not considered. Now significant input correlation only
exist between ’P1008’ and ‘P1012’, which is the same as we have observed in
steady-state cases in the previous section.

Figure 19: Posterior distribution and trace plots calibrated by 7T-FR,
covariance not considered Figure 20: Pair-wise joint distribution calibrated
by 7T-FR, covariance not considered Figure 21: Posterior distribution and
trace plots calibrated by 7T-FR, covariance is considered Figure 22: Pair-wise
joint distribution calibrated by 7T-FR, covariance is considered

In order to determine which result is more reasonable and should be used, in
Figure 23 we compare the model responses using different physical model
parameters: (1) $\theta_{prior}=1$ (black dotted line), (2) $\theta_{post1}$
when covariance is considered (blue line with circles), and (3)
$\theta_{post2}$ when covariance is not considered (orange dotted lines).
Experiment data (red circles) is also included for comparison.

Figure 23: Comparison of model responses with different calibration parameters

Comparing the results with and without covariance in Figure 23, we can see
that the result with covariance information tends to be closer to the
experiment during transition periods, when covariance is small (see Figure
18). While the result without covariance information simply finds the curve
that may lead to the smallest overall error level. For example in the left
figure, the curve without covariance chooses to “fit” the tail part of the
time series, but fails to follow the transition period when void fraction is
increasing at around 160s. So now we can better understand the role of the
covariance matrix in Bayesian calibration of time series data: it encourages
us to put more “weight” on the overall shape of the time series rather than
focusing on reducing averaged error.

###  5.3 Using the Hierarchical Framework for Time-Dependent Problems

We have demonstrated the effects of covariance information on the results of
calibration through one transient case. However, as we see in Figure 23, the
large model discrepancy may lead to significant over-fitting issue: the
posterior distribution of $\bm{\theta}$ obtained from one transient would not
help achieve better agreement between simulation and experiment for other
transient cases. In light of this situation, we propose to use the
hierarchical model to avoid over-fitting.

Similar to the hierarchical structure shown before, we assume that calibration
parameter $\bm{\theta}$ can be different for each group. In this scenario,
each group represents one transient case, and there are 9 transients in all,
we use 5 of them as the training (calibration) set and the other 4 of them as
the testing set. Five transients randomly selected for training are:
‘5T-FR’,‘5T-PI’,‘6T-FR’,‘6T-TI’,‘7T-FR’. The priors of parameters are defined
below:

$\displaystyle\mu\\_P1008$ $\displaystyle\sim\mathbf{Unif}(0,5)$
$\displaystyle\sigma\\_P1008$ $\displaystyle\sim\mathbf{Unif}(0,1)$
$\displaystyle P1008$
$\displaystyle\sim\mathcal{N}(\mu\\_P1008,\sigma\\_P1008)$

Using the NUTS method we are able to obtain converged posterior samples for
all the hyperparameters (shared parameter) and per-group parameters. Figure 24
shows the posterior distribution for all hyperparameters, and Figure 25 shows
the posterior for each per-group parameter and their trace plots. The cluster-
specific parameter have relatively large differences, as we can see in Figure
25. This is possibly due to the large model discrepancy in some cases, some
parameter are overfitting the model discrepancy instead of showing their own
distributions. The hierarchical model actually takes the overfitting issue
into consideration by allowing these parameters to be in a common
distribution, and we can take this distribution as our updated knowledge about
those parameters, which is a conservative measure and can avoid overfitting.
For example, for parameter ‘P1008’, we can use
$\mathcal{N}(\overline{{\mu\\_P1008}},\overline{{\sigma\\_P1008}})$ as the
posterior, where $\overline{{\mu\\_P1008}}$ and $\overline{{\sigma\\_P1008}})$
can be obtained by their posterior samples.

Figure 24: Posterior distributions of shared parameters in the hierarchical
model Figure 25: Posterior distributions of per-group parameters in the
hierarchical model Figure 26: Comparison of error distributions in the
training set  Figure 27: Comparison of error distributions in the validation
set

Now we can validate the posteriors using the validation set, the validation
method is the same as we have been using in previous sections. Since multiple
time series are more difficult to visualize, we directly compare the
distribution of errors (Exp. void fraction $-$ Simu. void fraction) using
calibrated values off the calibration parameter $\bm{\theta}$. The comparison
results for training set and for the validation set are shown in Figure 26 and
27, respectively. We can see that the error distribution is slightly improved
in the training set, and is significantly improved in the testing set. The
fact that the calibrated model performance in the validation set is better
than in the training set demonstrates the capability to avoid overfitting in
hierarchical models.

##  6 Summary

This work tackles the IUQ problem for time-dependent datasets in nuclear TH
applications. A key advancement of this study was the incorporation of
covariance information in posterior distribution calculations, demonstrating
how this inclusion better aligns the simulation results with experimental
data, especially during transition phases of the transients. This method
emphasizes capturing the time series’ overall shape rather than simply
minimizing average errors.

The central role of surrogate models in Bayesian calibration is underscored,
highlighting their utility in simplifying the computationally intensive
simulations. Artificial Neural Networks proved particularly effective for
handling high-dimensional outputs and for facilitating derivative calculations
within the MCMC algorithms, outperforming other models when applied to certain
TH applications.

The research also identified the risk of over-fitting with traditional IUQ
approaches in individual transients. To mitigate this, a hierarchical model
was utilized and validated, which showed improved error distributions for both
training and validation sets. This confirms the hierarchical model’s
superiority in avoiding over-fitting, thus presenting a robust solution for
IUQ in time-dependent nuclear TH simulations.

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G.S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., Zheng, X., 2015.  TensorFlow: Large-scale machine learning on heterogeneous systems.  URL: https://www.tensorflow.org/. software available from tensorflow.org.
  * Baccou et al. (2019) Baccou, J., Zhang, J., Fillion, P., Damblin, G., Petruzzi, A., Mendizábal, R., Reventós, F., Skorek, T., Couplet, M., Iooss, B., et al., 2019.  Development of good practice guidance for quantification of thermal-hydraulic code model input uncertainty.  Nuclear Engineering and Design 354, 110173.
  * Baccou et al. (2020) Baccou, J., Zhang, J., Fillion, P., Damblin, G., Petruzzi, A., Mendizábal, R., Reventós, F., Skorek, T., Couplet, M., Iooss, B., et al., 2020.  Sapium: A generic framework for a practical and transparent quantification of thermal-hydraulic code model input uncertainty.  Nuclear Science and Engineering 194, 721–736.
  * Borowiec et al. (2017) Borowiec, K., Wang, C., Kozlowski, T., Brooks, C.S., 2017.  Uncertainty quantification for steady-state psbt benchmark using surrogate models.  Transactions of the American Nuclear Society 117, 119–122.
  * Buitinck et al. (2013) Buitinck, L., Louppe, G., Blondel, M., Pedregosa, F., Mueller, A., Grisel, O., Niculae, V., Prettenhofer, P., Gramfort, A., Grobler, J., Layton, R., VanderPlas, J., Joly, A., Holt, B., Varoquaux, G., 2013.  API design for machine learning software: experiences from the scikit-learn project, in: ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pp. 108–122.
  * Chen et al. (2019) Chen, S., Kong, N., Sun, X., Meng, H., Li, M., 2019.  Claims data-driven modeling of hospital time-to-readmission risk with latent heterogeneity.  Health care management science 22, 156–179.
  * Chen et al. (2018) Chen, S., Lu, L., Xiang, Y., Lu, Q., Li, M., 2018.  A data heterogeneity modeling and quantification approach for field pre-assessment of chloride-induced corrosion in aging infrastructures.  Reliability Engineering & System Safety 171, 123–135.
  * Chen et al. (2023) Chen, S., Wu, J., Hovakimyan, N., Yao, H., 2023.  Recontab: Regularized contrastive representation learning for tabular data.  arXiv preprint arXiv:2310.18541 .
  * Dong et al. (2021a) Dong, G., Cai, L., Datta, D., Kumar, S., Barnes, L.E., Boukhechba, M., 2021a.  Influenza-like symptom recognition using mobile sensing and graph neural networks, in: Proceedings of the conference on health, inference, and learning, pp. 291–300.
  * Dong et al. (2021b) Dong, G., Tang, M., Cai, L., Barnes, L.E., Boukhechba, M., 2021b.  Semi-supervised graph instance transformer for mental health inference, in: 2021 20th IEEE International Conference on Machine Learning and Applications (ICMLA), IEEE. pp. 1221–1228.
  * Gal and Ghahramani (2016) Gal, Y., Ghahramani, Z., 2016.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in: international conference on machine learning, pp. 1050–1059.
  * Ghione et al. (2023) Ghione, A., Sargentini, L., Damblin, G., Fillion, P., Baccou, J., Sueur, R., Iooss, B., Petruzzi, A., Zeng, K., Zhang, J., Havet, M., Mendizábal, R., Skorek, T., Wu, X., Freixa Terradas, J., Adorni, M., 2023.  Applying the sapium guideline for input uncertainty quantification: the atrium project, in: Proceedings of the 20th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-20), Washington, D.C., USA, August 20-25, 2023.
  * Hecht-Nielsen (1992) Hecht-Nielsen, R., 1992.  Theory of the backpropagation neural network, in: Neural networks for perception. Elsevier, pp. 65–93.
  * Hu et al. (2022) Hu, J., Xu, Y., Tang, Z., 2022.  Dan-pd: Domain adaptive network with parallel decoder for polyp segmentation.  Computerized Medical Imaging and Graphics 101, 102124.
  * Kennedy and O’Hagan (2001) Kennedy, M.C., O’Hagan, A., 2001.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63, 425–464.
  * Li et al. (2023a) Li, Z., Kong, X., Zhang, Y., 2023a.  Exploring factors associated with crossing assertiveness of pedestrians at unsignalized intersections.  Transportation Research Record 2677, 182–198.
  * Li et al. (2023b) Li, Z., Tang, Z., Hu, J., Wang, X., Jia, D., Zhang, Y., 2023b.  Nst: A nuclei segmentation method based on transformer for gastrointestinal cancer pathological images.  Biomedical Signal Processing and Control 84, 104785.
  * Li et al. (2023c) Li, Z., Wei, Z., Zhang, Y., Kong, X., Ma, C., 2023c.  Applying an interpretable machine learning framework to study mobility inequity in the recovery phase of covid-19 pandemic.  Travel Behaviour and Society 33, 100621.
  * Liu and Bao (2023) Liu, Y., Bao, Y., 2023.  Intelligent monitoring of spatially-distributed cracks using distributed fiber optic sensors assisted by deep learning.  Measurement 220, 113418.
  * Liu et al. (2018) Liu, Y., Dinh, N., Sato, Y., Niceno, B., 2018.  Data-driven modeling for boiling heat transfer: using deep neural networks and high-fidelity simulation results.  Applied Thermal Engineering 144, 305–320.
  * Liu et al. (2021) Liu, Y., Liu, L., Yang, L., Hao, L., Bao, Y., 2021.  Measuring distance using ultra-wideband radio technology enhanced by extreme gradient boosting decision tree (xgboost).  Automation in Construction 126, 103678.
  * Ma et al. (2022) Ma, C., Peng, Y., Wu, L., Guo, X., Wang, X., Kong, X., 2022.  Application of machine learning techniques to predict the occurrence of distraction-affected crashes with phone-use data.  Transportation research record 2676, 692–705.
  * Meng et al. (2022) Meng, Y., Wu, L., Ma, C., Guo, X., Wang, X., 2022.  A comparative analysis of intersection hotspot identification: Fixed vs. varying dispersion parameters in negative binomial models.  Journal of Transportation Safety & Security 14, 305–322.
  * Perez et al. (2011) Perez, M., Reventos, F., Batet, L., Guba, A., Tóth, I., Mieusset, T., Bazin, P., De Crécy, A., Borisov, S., Skorek, T., et al., 2011.  Uncertainty and sensitivity analysis of a lbloca in a pwr nuclear power plant: Results of the phase v of the bemuse programme.  Nuclear Engineering and Design 241, 4206–4222.
  * Ramachandran et al. (2017) Ramachandran, P., Zoph, B., Le, Q.V., 2017.  Searching for activation functions.  arXiv preprint arXiv:1710.05941 .
  * Rasmussen (2004) Rasmussen, C.E., 2004.  Gaussian processes in machine learning, in: Advanced lectures on machine learning. Springer, pp. 63–71.
  * Reed and MarksII (1999) Reed, R., MarksII, R.J., 1999.  Neural smithing: supervised learning in feedforward artificial neural networks.  Mit Press.
  * Rubin et al. (2012) Rubin, A., Schoedel, A., Avramova, M., Utsuno, H., Bajorek, S., Velazquez-Lozada, A., 2012.  OECD/NRC Benchmark Based on NUPEC PWR Sub-channel and Bundle Test (PSBT). Volume I: Experimental Database and Final Problem Specifications.  Technical Report. Organisation for Economic Co-Operation and Development.
  * Santner et al. (2013) Santner, T.J., Williams, B.J., Notz, W.I., 2013.  The design and analysis of computer experiments.  Springer Science & Business Media.
  * Shlens (2014) Shlens, J., 2014.  A tutorial on principal component analysis.  arXiv preprint arXiv:1404.1100 .
  * Skorek et al. (2019) Skorek, T., de Crécy, A., Kovtonyuk, A., Petruzzi, A., Mendizábal, R., de Alfonso, E., Reventós, F., Freixa, J., Sarrette, C., Kyncl, M., et al., 2019.  Quantification of the uncertainty of the physical models in the system thermal-hydraulic codes–premium benchmark.  Nuclear Engineering and Design 354, 110199.
  * Snoek et al. (2012) Snoek, J., Larochelle, H., Adams, R.P., 2012.  Practical bayesian optimization of machine learning algorithms, in: Advances in neural information processing systems, pp. 2951–2959.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R., 2014.  Dropout: a simple way to prevent neural networks from overfitting.  The journal of machine learning research 15, 1929–1958.
  * Stewart et al. (2008) Stewart, L.M., Dance, S., Nichols, N., 2008.  Correlated observation errors in data assimilation.  International journal for numerical methods in fluids 56, 1521–1527.
  * Stewart et al. (2013) Stewart, L.M., Dance, S.L., Nichols, N.K., 2013.  Data assimilation with correlated observation errors: experiments with a 1-d shallow water model.  Tellus A: Dynamic Meteorology and Oceanography 65, 19546.
  * Stripling et al. (2013) Stripling, H., McClarren, R., Kuranz, C., Grosskopf, M., Rutter, E., Torralva, B., 2013.  A calibration and data assimilation method using the bayesian mars emulator.  Annals of Nuclear Energy 52, 103–112.
  * Tiedeman and Green (2013) Tiedeman, C.R., Green, C.T., 2013.  Effect of correlated observation error on parameters, predictions, and uncertainty.  Water Resources Research 49, 6339–6355.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Jarrod Millman, K., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C., Polat, İ., Feng, Y., Moore, E.W., Vand erPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., Contributors, S…, 2020.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods 17, 261–272.  doi:https://doi.org/10.1038/s41592-019-0686-2.
  * Waller et al. (2016) Waller, J.A., Dance, S.L., Nichols, N.K., 2016.  Theoretical insight into diagnosing observation error correlations using observation-minus-background and observation-minus-analysis statistics.  Quarterly Journal of the Royal Meteorological Society 142, 418–431.
  * Wang (2020) Wang, C., 2020.  A hierarchical Bayesian calibration framework for quantifying input uncertainties in thermal-hydraulics simulation models.  Ph.D. thesis. University of Illinois at Urbana-Champaign.
  * Wang et al. (2018a) Wang, C., Wu, X., Borowiec, K., Kozlowski, T., 2018a.  Bayesian calibration and uncertainty quantification for trace based on psbt benchmark.  Transactions of the American Nuclear Society 118, 419–422.
  * Wang et al. (2017a) Wang, C., Wu, X., Kozlowski, T., 2017a.  Sensitivity and uncertainty analysis of trace physical model parameters based on psbt benchmark using gaussian process emulator.  Proc. 17th Int. Topl. Mtg. Nuclear Reactor Thermal Hydraulics (NURETH-17) , 3–8.
  * Wang et al. (2017b) Wang, C., Wu, X., Kozlowski, T., 2017b.  Surrogate-based inverse uncertainty quantification of trace physical model parameters using steady-state psbt void fraction data.  Proc. 17th Int. Topl. Mtg. Nuclear Reactor Thermal Hydraulics (NURETH-17) , 3–8.
  * Wang et al. (2018b) Wang, C., Wu, X., Kozlowski, T., 2018b.  Surrogate-based bayesian calibration of thermal-hydraulics models based on psbt time-dependent benchmark data, in: Proceedings of the ANS Best Estimate Plus Uncertainty International Conference (BEPU-2018), Real Collegio, Lucca, Italy, May 13-19, 2018.
  * Wang et al. (2019a) Wang, C., Wu, X., Kozlowski, T., 2019a.  Gaussian process–based inverse uncertainty quantification for trace physical model parameters using steady-state psbt benchmark.  Nuclear Science and Engineering 193, 100–114.
  * Wang et al. (2019b) Wang, C., Wu, X., Kozlowski, T., 2019b.  Inverse uncertainty quantification by hierarchical bayesian inference for trace physical model parameters based on bfbt benchmark, in: Proc. 18th Int. Topl. Mtg. Nuclear Reactor Thermal Hydraulics (NURETH-18), Portland, OR, USA.
  * Wang et al. (2023a) Wang, C., Wu, X., Kozlowski, T., 2023a.  Inverse uncertainty quantification by hierarchical bayesian modeling and application in nuclear system thermal-hydraulics codes.  arXiv preprint arXiv:2305.16622 .
  * Wang et al. (2023b) Wang, C., Wu, X., Xie, Z., Kozlowski, T., 2023b.  Scalable inverse uncertainty quantification by hierarchical bayesian modeling and variational inference.  Energies 16, 7664.
  * Wang et al. (2023c) Wang, H., Fu, T., Du, Y., Gao, W., Huang, K., Liu, Z., Chandak, P., Liu, S., Van Katwyk, P., Deac, A., et al., 2023c.  Scientific discovery in the age of artificial intelligence.  Nature 620, 47–60.
  * Wu et al. (2022) Wu, J., Tao, R., Zhao, P., Martin, N.F., Hovakimyan, N., 2022.  Optimizing nitrogen management with deep reinforcement learning and crop simulations, in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 1712–1720.
  * Wu and Kozlowski (2017) Wu, X., Kozlowski, T., 2017.  Inverse uncertainty quantification of reactor simulations under the bayesian framework using surrogate models constructed by polynomial chaos expansion.  Nuclear Engineering and Design 313, 29–52.
  * Wu et al. (2018) Wu, X., Kozlowski, T., Meidani, H., 2018.  Kriging-based inverse uncertainty quantification of nuclear fuel performance code bison fission gas release model using time series measurement data.  Reliability Engineering & System Safety 169, 422–436.
  * Wu et al. (2017a) Wu, X., Mui, T., Hu, G., Meidani, H., Kozlowski, T., 2017a.  Inverse uncertainty quantification of trace physical model parameters using sparse gird stochastic collocation surrogate model.  Nuclear Engineering and Design 319, 185–200.
  * Wu et al. (2017b) Wu, X., Wang, C., Kozlowski, T., 2017b.  Global sensitivity analysis of trace physical model parameters based on bfbt benchmark, in: Proceedings of the MC-2017, International Conference on Mathematics Computational Methods Applied to Nuclear Science Engineering, Jeju, Korean.
  * Wu et al. (2017c) Wu, X., Wang, C., Kozlowski, T., 2017c.  Kriging-based surrogate models for uncertainty quantification and sensitivity analysis, in: Proceedings of the MC-2017, International Conference on Mathematics Computational Methods Applied to Nuclear Science Engineering, Jeju, Korean.
  * Xue et al. (2022) Xue, B., Jiao, Y., Kannampallil, T., Fritz, B., King, C., Abraham, J., Avidan, M., Lu, C., 2022.  Perioperative predictions with interpretable latent representation, in: Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 4268–4278.
  * Xue et al. (2021) Xue, B., Li, D., Lu, C., King, C.R., Wildes, T., Avidan, M.S., Kannampallil, T., Abraham, J., 2021.  Use of machine learning to develop and evaluate models using preoperative and intraoperative data to identify risks of postoperative complications.  JAMA network open 4, e212240–e212240.
  * Yurko et al. (2015) Yurko, J.P., Buongiorno, J., Youngblood, R., 2015.  Demonstration of emulator-based bayesian calibration of safety analysis codes: theory and formulation.  Science and Technology of Nuclear Installations 2015.


# A short proof of the existence of a minor-universal countable planar graph

George Kontogeorgiou

An element of a class of graphs is _minor-universal_ for the class if every
other element of the class is its minor. In [2], Diestel and Kühn constructed
a minor-universal countable planar graph. Their construction is rather
involved. The purpose of this note is to present a significantly simpler
construction and proof, based on a suggestion of Georgakopoulos ([3], p.7) to
modify the construction of [4]. The idea is as follows:

  1. i.

$\mathcal{G}_{0}$ is a copy of $C_{3}$ embedded in $\mathbb{S}^{2}$;

  2. ii.

for $n>0$, $\mathcal{G}_{n}$ is the plane graph obtained from
$\mathcal{G}_{n-1}$ in the following way: in each face $F$ of
$\mathcal{G}_{n-1}$ we place a vertex $v_{F}$, we connect $v_{F}$ to each
vertex $u\in\partial F$ with an edge, and we subdivide each edge
$\\{v_{F},u\\}$ exactly $n$ times;

  3. iii.

we define $\mathcal{G}:=\bigcup_{n=0}^{\infty}\mathcal{G}_{n}$.

Figure 1: $\mathcal{G}_{2}$ is formed by gluing two copies of this graph at
the outer triangle.

In this note I prove:

######  Theorem 1.

The graph $\mathcal{G}$ is minor-universal for countable planar graphs.

By its construction, $\mathcal{G}$ is countable and planar. It remains to show
that it has every other countable planar graph as a minor. This follows
directly from the following two lemmata:

######  Lemma 1.

Every countable planar graph is a minor of a countable, sub-cubic, 2-connected
planar graph.

###### Proof.

Let $G\in\mathcal{C}$. It is simple to construct a countable planar graph
$G^{\prime}\succ G$ in $\mathcal{C}$ that is 2-connected and has finite
maximum degree by blowing up each vertex to a locally finite tree and taking
the fattening, for details see ([3], Chapter 4). We blow up each vertex of
$G^{\prime}$ to a finite trivalent tree to get the promised graph
$G^{\prime\prime}$. ∎

######  Lemma 2.

Every countable, sub-cubic, 2-connected planar graph is a minor of
$\mathcal{G}$.

The rest of this note is devoted to proving Lemma 2.

## Proof of Lemma 2

An _ear decomposition_ of a countable graph $G$ is a sequence
$(G_{n})_{n\in\mathbb{N}}$ of subgraphs of $G$ such that:

  1. i.

$G_{0}$ is a cycle;

  2. ii.

for $n>0$, either $G_{n}=G_{n-1}$ or $G_{n}$ is a simple graph obtained from
$G_{n-1}$ by connecting two distinct vertices of $G_{n-1}$ with a new path
(called an _ear_ in this context);

  3. iii.

$G=\bigcup_{n=0}^{\infty}G_{n}$.

We first prove two claims that are used in our proof.

######  Claim 1.

Every 2-connected countable graph $G$ with at least three vertices admits an
ear decomposition.

###### Proof.

We order the edges of $G$. By the infinite version of Menger’s Theorem ([1]),
there exists a cycle $C$ that contains $e_{0}$. Set $G_{0}:=C$. Given
$G_{n-1}$, we determine $G_{n}$ as follows. If $G_{n-1}=G$, then set
$G_{n}:=G$. Otherwise, let $e$ be the edge outside of $G_{n-1}$ with the
smallest index. It is a well-known exercise that there exists a cycle $C_{e}$
in $G$ that contains both $e$ and $e_{0}$. Let $P_{e}$ be the subpath of
$C_{e}$ that meets $G_{n-1}$ only at its endpoints and contains $e$. Set
$G_{n}:=G_{n-1}\cup P_{e}$. To see that $G=\bigcup_{n=0}^{\infty}G_{n}$, note
that $e_{n}\in G_{n}$.

∎

######  Remark 1.

If $G$ is a plane graph, then the $n^{th}$ ear of its ear decomposition lies
in a face of $G_{n-1}$.

An $n$-diameter of $\mathcal{G}$ is a maximal path contained in
$\mathcal{G}_{n+1}\setminus\mathcal{G}_{n}$. An $n$-slice of $\mathcal{G}$ is
a subgraph that is bounded by (but does not contain the edges of) the boundary
of a face of $\mathcal{G}_{n}$ and perhaps one of the diameters of that face,
and does not contain $\mathcal{G}_{0}$. An $n$-piece is a subgraph of
$\mathcal{G}$ of the form $\overline{F}\cap\mathcal{G}$ or
$(\overline{F}\cap\mathcal{G})\setminus\mathcal{G}_{n}$ for some face $F$ of
$\mathcal{G}_{n}$. Recall that, given a separator $S$ of a graph $G$ and a
connected component $K\subseteq G\setminus S$, the _torso_ of $K$ is the
induced subgraph of $G$ on the vertices $V(K)\cup S$.

Figure 2: A purple $2$-diameter bounding a grey area that contains a
$2$-slice.

######  Claim 2.

Every $n$-slice is $(n+3)$-connected.

###### Proof.

Let $J$ be an $n$-slice and $S$ be a separator of $J$. We say that a connected
component is $N$-_coarse_ for $N\geq n$ if its torso in $J$ is the union of a
finite set of $N$-pieces. Since $S$ is finite, it is entirely contained in
some $\mathcal{G}_{N}$ for $N\geq n$, so it splits $J$ into $N$-coarse
components. Therefore $S$ must have order at least $N+3\geq n+3$, since $N+3$
is half the length of an $N$-diameter rounded up, which is the shortest
boundary that an $N$-coarse component can have with its complement in $J$. ∎

Recall that an inflated copy (briefly denoted as $IG$) of a graph $G=(V,E)$ is
a graph $H$ together with a map $\phi:G\rightarrow H$ such that for every
$v\in V$ and $e\in E$ with ends $\\{x,y\\}$, the following are true:

  1. i.

$\phi(v)$ is a connected subgraph of $H$;

  2. ii.

$\phi(e)$ is an open path between $\phi(x)$ and $\phi(y)$;

  3. iii.

the images of the elements of $G$ are pairwise disjoint;

  4. iv.

$\phi(G)=H$.

A graph $G$ is a minor of a graph $G^{\prime}$ if and only if $G^{\prime}$
contains an $IG$ as a subgraph. We are now ready to prove Lemma 2:

###### Proof.

Let $G$ be a plane graph as in the statement of Lemma 2. We will construct a
subgraph $H$ of $\mathcal{G}$ that is an $IG$. In particular, we will define
$H$ as the union of a sequence $H_{n}$ of $IG_{n}$, where
$(G_{n})_{n\in\mathbb{N}}$ is an ear decomposition of $G$ (Claim 1), such that
for every $n\in\mathbb{N}$, $v\in V(G)$, and $e\in E(G)$, we have for the
corresponding maps that $\phi_{n}(v)\subseteq\phi_{n+1}(v)$ and
$\phi_{n}(e)=\phi_{n+1}(e)$. This will allow us to define $\phi:G\rightarrow
H$ by setting $\phi:=\bigcup_{n=0}^{\infty}\phi_{n}$; the reader is invited to
check that this yields an IG. To satisfy these constraints, we also require a
consistent way to correspond abstract faces of $G_{n}$ to slices of
$\mathcal{G}$. To keep notation simple, we overload our maps $\phi_{n}$ to
also perform this function.

The following definition is key in our analysis: given a face $F$ of $G_{n}$,
we denote by $V_{F}$ the set of vertices of $\partial F$ that are ends of ears
contained in $F$. Note that, since $G$ is sub-cubic, each vertex can be the
end of at most one ear.

We construct inductively the sequence $\phi_{n}$ with the additional useful
requirements that for each $v\in V(G_{n})$, $\phi_{n}(v)$ is a path, and that
for each face $F$ of $G_{n}$, the slice $\phi_{n}(F)$ is $|V_{F}|$-connected
and intersects $H_{n}$ only at one end vertex of each of the paths
$\\{\phi_{n}(v)|v\in\partial F\\}$ and in the correct cyclic order.

For the base case, suppose that $G_{0}$ has length $l_{0}\geq 3$. We consider
two disjoint face cycles $C_{0}$ and $C^{\prime}_{0}$ of the graph
$\mathcal{G}_{l_{0}}$, contained in the same ($l_{0}-2$)-piece, noting that
each has length $2l_{0}+3>l_{0}$. By Menger’s Theorem and Claim 2, we may also
consider $l_{0}$ disjoint paths inside the piece joining $C_{0}$ and
$C^{\prime}_{0}$. We define $\phi_{0}$ by mapping the vertices of $G_{0}$ to
these paths in the correct cyclic order, each edge of $G_{0}$ to an
appropriate (open) subpath of $C_{0}$, and the two faces $F_{0}$ and
$F^{\prime}_{0}$ of $G_{0}$ to the two $l_{0}$-slices $J_{0}$ and
$J^{\prime}_{0}$ bounded by $C_{0}$ and $C^{\prime}_{0}$, respectively. Both
$J_{0}$ and $J^{\prime}_{0}$ are $(l_{0}+3)$-connected, and
$l_{0}+3>l_{0}=|V_{F_{0}}|+|V_{F^{\prime}_{0}}|$, so $F_{0}$ is mapped to a
$|V_{F_{0}}|$-connected slice and $F^{\prime}_{0}$ is mapped to a
$|V_{F^{\prime}_{0}}|$-connected slice.

For the inductive step, suppose that we have defined $\phi_{n-1}$ so that each
vertex $v$ is mapped to a path $\phi_{n-1}(v)$ and each face $F$ is mapped to
a $|V_{F}|$-connected slice $\phi_{n-1}(F)$ that intersects $H_{n-1}$ only at
one end vertex of each of the paths $\\{\phi_{n-1}(v)|v\in\partial F\\}$ and
in the correct cyclic order.

Let $F$, of boundary length say $l$, be the unique (Remark 1) face of
$G_{n-1}$ such that the $n^{th}$ ear $G_{n}\setminus G_{n-1}$, of length say
$L$, lies within $F$, splitting it to two faces $F^{\prime}$ and
$F^{\prime\prime}$. Let $f$ be a face of $\mathcal{G}_{N}$ for some $N>n$ such
that $f$ is bounded by the slice $\phi_{n-1}(F)$, has diameter length at least
$L$, and all its slices have connectivity at least $2L+l$ (e.g. take
$N=2L+l-3$) (Claim 2).

Figure 3: A step in the construction of $\phi$.

By the inductive hypothesis, $\phi_{n-1}(F)$ is $|V_{F}|$-connected. By
Menger’s Theorem, there exist disjoint paths $(P_{v})_{v\in V_{F}}$ in
$\phi_{n-1}(F)$ such that $P_{v}$ begins at the vertex
$\phi_{n-1}(v)\cap\phi_{n-1}(F)$ and ends at some $p(v)\in\partial f$. In
particular, since $\mathcal{G}$ is a plane graph, $p$ preserves the cyclic
order of $\phi_{n-1}(V_{F})$. We also denote by $D$ the diameter of $f$ with
ends $p(x)$ and $p(y)$, where $x$ and $y$ are the ends of the ear
$G_{n}\setminus G_{n-1}$.

For $v\in V_{F}$ we set $\phi_{n}(v)=\phi_{n-1}(v)\cup P_{v}$. The vertices of
the ear $G_{n}\setminus G_{n-1}$ we map arbitrarily (but in the correct order)
to single vertices of $D$. Each edge of the ear is mapped to the subpath of
$D$ between the images of its ends. The faces $F^{\prime}$ and $F$ are mapped
to the corresponding slices induced on $f$ by $D$. Every other vertex, edge,
or face of $G_{n}$ has the same image under $\phi_{n}$ as under $\phi_{n-1}$.

Firstly, it follows from the construction that $\phi_{n}:G_{n}\rightarrow
H_{n}$ is an $IG_{n}$.

Moreover, for $v\in V_{F}$, by the inductive hypothesis, $\phi_{n-1}(v)$ is a
path. Also, $P_{v}$ is a path and, again by the inductive hypothesis,
$\phi_{n-1}(v)\cap P_{v}$ is a single vertex. Therefore, $\phi_{n}(v)$ is a
path. Trivially, the same holds for vertices in $V(G_{n})\setminus V_{F}$.

Additionally, both $F^{\prime}$ and $F^{\prime\prime}$ are mapped to slices of
connectivity $2L+l$. Since
$|V_{F^{\prime}}|+|V_{F^{\prime\prime}}|\leq|\partial F^{\prime}|+|\partial
F^{\prime\prime}|=2L+l$, both of these slices are adequately connected. We
have also ensured that they intersect the appropriate vertex images under
$\phi_{n}$ only at one end vertex each and in the correct cyclic order. These
same properties hold, by the inductive hypothesis, for every other face image
of $G_{n}$. This concludes the induction, and the proof of Lemma 2.

∎

######  Remark 2.

The above construction is quite manageable, and perhaps can be modified to
tackle the problem of minor-universal elements for 2-complexes embeddable in
$\mathbb{S}^{3}$, mentioned in ([3], p.12).

## References

  * [1] Halin, R. (1974). A note on Menger’s theorem for infinite locally finite graphs. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg. 40:111-114.
  * [2] Diestel, R. and Kühn, D. (1999). A universal planar graph under the minor relation. Journal of Graph Theory. 32:191-206.
  * [3] Georgakopoulos, A. (2022). On graph classes with minor-universal elements. Preprint.
  * [4] Georgakopoulos, A. (2006). Infinite highly connected planar graphs of large girth. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg. 76(1):235-245.