# Propagating uncertainty in a network of energy models ††thanks: This work
was funded by the EPSRC National Centre for Energy Systems Integration through
the flex fund award FFC3-008, and by the Alan Turing Institute grant ‘Managing
Uncertainty in Government Modelling’. VV, NS and HPW were contracted by the
Turing Institute when this work was carried out, and JQS, PGC and CJD are
Turing Fellows at the Institute.

Victoria Volodina  Clinical Operational
Research Unit
University College London
London, UK
<EMAIL_ADDRESS>Nikki Sonenberg  Heilbronn Institute of
Mathematical Research
University of Bristol
Bristol, UK
<EMAIL_ADDRESS>Jim Q. Smith
Peter G. Challenor  Department of Statistics
Warwick University
Coventry, UK
<EMAIL_ADDRESS>Department of Mathematics
University of Exeter
Exeter, UK
<EMAIL_ADDRESS>Chris J. Dent  School of Mathematics
University of Edinburgh
Edinburgh, UK
<EMAIL_ADDRESS>Henry P. Wynn  Department of Statistics
London School of Economics
London, UK
<EMAIL_ADDRESS>
###### Abstract

Computer models are widely used in decision support for energy systems
operation, planning and policy. A system of models is often employed, where
model inputs themselves arise from other computer models, with each model
being developed by different teams of experts. Gaussian Process emulators can
be used to approximate the behaviour of complex, computationally intensive
models and used to generate predictions together with a measure of uncertainty
about the predicted model output. This paper presents a computationally
efficient framework for propagating uncertainty within a network of models
with high-dimensional outputs used for energy planning. We present a case
study from a UK county council considering low carbon technologies to
transform its infrastructure to reach a net-zero carbon target. The system
model considered for this case study is simple, however the framework can be
applied to larger networks of more complex models.

###### Index Terms:

energy systems, decision support, surrogate, Gaussian processes, uncertainty
propagation

##  I Introduction

Computer models are widely used in decision support for energy systems
operation, planning and policy. A network of models is often employed, where
model inputs themselves arise from other computer models, with each model
being developed by different experts from across various disciplines. To
perform inferences on a single computer model such as calibration, prediction,
uncertainty and sensitivity analysis, we require multiple simulation runs.
However, as the computer model can be very expensive to operate, we need to
acknowledge the uncertainty about its outputs at unseen input parameter
values, namely code uncertainty [1]. This is relevant for the analysis of
networks of models, since the outputs together with uncertainties from the
first layer computer models determine the uncertain inputs in the second layer
computer models and so on. For principled decision support under UK government
guidelines [2], the uncertainties associated with individual models’ outputs
need to be explicitly quantified and propagated across the network. Other
related frameworks for the coupling of probabilistic systems include [3].

Gaussian Process (GP) emulators are commonly used to approximate the behaviour
of complex, computationally intensive computer models; an emulator provides
predictions and quantifies uncertainty in the state of knowledge regarding the
behaviour of the model including that arising from only a limited number of
runs being possible. GP emulators have been widely used as _surrogates_ to
complex computer models in climate and environmental studies [4, 5, 6] and
electricity prices [7].

This paper presents a framework for propagating uncertainty within a network
of models with high-dimensional outputs (time-series) using recent
methodological advances in the theory of networks of GP emulators [8, 9, 10].
First, we employ principal component analysis (PCA) to reduce the dimension of
the output space by projecting the high-dimensional output onto a low-
dimensional basis, we then specify a univariate GP emulator for the individual
coefficients of a low-dimensional basis. The efficient propagation of
uncertainty between models is achieved by passing the first and second moments
between the probability models for the basis vectors’ coefficients obtained
and projecting back to reconstruct the mean and variance for the original
model output.

This methodology is applied to a case study for a UK county council
considering low carbon technology options to transform the infrastructure at a
facility in order to reach a carbon zero target in 2050 [11]. Under the Public
Decarbonisation Scheme [12], the council plans to replace the gas boiler with
an electric-powered ground source heat pump (GSHP) to supply heating. The
simple network model is used to demonstrate how projections of operational
costs of the gas boiler can be generated, however the framework is scalable
and can be applied to larger networks of more complex models.

In Section II we describe the GP methodology, and in Section III we describe
the individual computer models used in energy systems planning. We present
results for the case study in Section IV and in Section V we provide a
discussion.

##  II Methodology

In Sections II-A & II-B we define the univariate GP model and multivariate GP
models and use these to assemble the linked GP emulator in Section II-C.

###  II-A Gaussian Process model

Let $\boldsymbol{\mathrm{x}}=(x_{1},\dots,x_{p})\in\mathbb{R}^{p}$ be a
$p$-dimensional vector of inputs and $f(\boldsymbol{\mathrm{x}})$ be the
scalar-valued output that represents the process of interest. The GPs are
fully specified by mean function $\mu(\cdot)$ and covariance function
$\sigma^{2}r(\cdot,\cdot;\boldsymbol{\delta})$, where $\boldsymbol{\delta}$
and $\sigma^{2}$ correspond to the vector of correlation length parameters and
variance parameter respectively.

Suppose we observe $n$ realisations
$F^{\mathcal{T}}=(f(\boldsymbol{\mathrm{x}}_{1}),\dots,f(\boldsymbol{\mathrm{x}}_{n}))$
at design points
$X^{\mathcal{T}}=(\boldsymbol{\mathrm{x}}_{1},\dots,\boldsymbol{\mathrm{x}}_{n})$.
The emulator is fitted based on an ensemble of runs of the model $f$, denoted
by $\mathcal{D}=\\{X^{\mathcal{T}},F^{\mathcal{T}}\\}$, using a Bayesian
approach with a non-informative prior for parameters $\sigma^{2}$ and
$\boldsymbol{\delta}$ [13].

###  II-B Multivariate Gaussian Process model

For computer models with multivariate outputs given by an $l$-dimensional
vector
$\boldsymbol{f}(\boldsymbol{\mathrm{x}})=(f_{1}(\boldsymbol{\mathrm{x}}),\dots,f_{l}(\boldsymbol{\mathrm{x}}))$,
for example, time series output, the data can be projected onto a low-
dimensional basis using principal components [14], and then independent GP
models are specified for coefficients of this basis.

The $n$ computer model simulations are stored in an $l\times n$ matrix
$\boldsymbol{F}=(\boldsymbol{f}(\boldsymbol{\mathrm{x}}_{1}),\dots,\boldsymbol{f}(\boldsymbol{\mathrm{x}}_{n}))$.
After subtracting out the mean simulation, $\boldsymbol{\mu}$, and scaling the
simulation output to obtain the centred ensemble, $\boldsymbol{F}_{\mu}$, the
singular value decomposition for the centred ensemble can be written as

$\boldsymbol{F}_{\mu}^{T}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{T}.$
(1)

We then obtain the principal component basis, denoted by
$\boldsymbol{\Gamma}=(\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{n-1})$,
which are the first $n-1$ columns of $\boldsymbol{V}$ and where each
individual basis vector $\boldsymbol{\gamma}_{i}$ has length $l$.

Given this basis, $\boldsymbol{f}(\cdot)$ can be written as a linear
combination of the basis vectors:

$\displaystyle\boldsymbol{f}(\boldsymbol{\mathrm{x}})-\boldsymbol{\mu}$
$\displaystyle=\sum_{i=1}^{n-1}\boldsymbol{\gamma}_{i}c_{i}(\boldsymbol{\mathrm{x}})+\boldsymbol{\epsilon},$
(2)
$\displaystyle=\boldsymbol{\Gamma}\boldsymbol{c}(\boldsymbol{\mathrm{x}})+\boldsymbol{\epsilon},$
(3)

where $c_{i}(\boldsymbol{\mathrm{x}})$ is the coefficient for basis vector
$\boldsymbol{\gamma}_{i}$ and $\boldsymbol{\epsilon}$ is a residual vector.

Exploiting the fact that the basis vectors of $\boldsymbol{\Gamma}$ are
orthogonal, we can fit univariate Gaussian process emulators for the
coefficients $c_{i}(\cdot)$ for each basis vector separately. We choose the
first $q$ vectors and denote the truncated basis by
$\boldsymbol{\Gamma}_{q}=(\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{q})$.
Conditioned on the projected ensemble, we emulate the coefficients for the $q$
basis vectors and obtain the posterior distribution for the reconstruction of
$\boldsymbol{f}(\cdot)$.

###  II-C Linked Gaussian Process emulators

To perform the coupling of a network of feed-forward computer models, we
present the linked GP emulator methodology used for uncertainty
quantification. Consider a simple system, where a $d$-dimensional output
produced by computer models in the first layer feed into a computer model in
the second layer [8].

For $i=1,\dots,d$ define a GP emulator to approximate the outputs produced by
individual computer models

$\displaystyle
W_{i}|\boldsymbol{\mathrm{x}}_{i}\sim\text{GP}(\mu_{i}(\cdot),\sigma_{i}^{2}r(\cdot,\cdot;\boldsymbol{\delta}_{i})),$
(4)

and

$\displaystyle
Y|W_{1},\dots,W_{d},\boldsymbol{z}\sim\text{GP}(\mu_{y}(\cdot),\sigma_{y}^{2}r(\cdot,\cdot;\boldsymbol{\delta}_{y})),$
(5)

where $\boldsymbol{\mathrm{x}}_{1},\dots,\boldsymbol{\mathrm{x}}_{d}$ and
$\boldsymbol{z}$ are vectors of input parameters. The outputs of the computer
models are treated as random variables, and the relationships between these
are shown in Figure 1.

$\boldsymbol{\mathrm{x}}_{1}$$\boldsymbol{\mathrm{x}}_{2}$$W_{1}$$W_{2}$$Y$$\boldsymbol{\mathrm{x}}_{d}$$W_{d}$$\boldsymbol{z}$$\dots$
Figure 1: A feed-forward graph depicting the relationship between random
variables $W_{1},W_{2},\dots,W_{d}$ and $Y$.

The distribution of
$Y|\boldsymbol{\mathrm{x}}_{1},\dots,\boldsymbol{\mathrm{x}}_{d},\boldsymbol{z}$
can be written as,

$\displaystyle
p(y|\boldsymbol{\mathrm{x}}_{1},..,\boldsymbol{\mathrm{x}}_{d},\boldsymbol{z})=\int
p(y|\boldsymbol{w},\boldsymbol{z})p(\boldsymbol{w}|\boldsymbol{\mathrm{x}}_{1},..,\boldsymbol{\mathrm{x}}_{d})d\boldsymbol{w},$

where $\boldsymbol{w}=(w_{1},\dots,w_{d})^{T}$. However,
$(y|\boldsymbol{\mathrm{x}}_{1},\dots,\boldsymbol{\mathrm{x}}_{d},\boldsymbol{z})$
is neither analytically tractable nor Gaussian in general.

The first two moments can be computed by Monte Carlo samples under the
assumption that the densities inside the integral are Gaussian [9]. Further,
under some mild conditions, the first two moments can be calculated
analytically [8, 10]. In particular, suppose we observe realisations of $y$
and $\boldsymbol{w}$, denoted by $Y^{\mathcal{T}}$ and $W^{\mathcal{T}}$, at
design (training) sets $X^{\mathcal{T}}$ and $Z^{\mathcal{T}}$, let
$\mathcal{D}=\\{X^{\mathcal{T}},Z^{\mathcal{T}},W^{\mathcal{T}},Y^{\mathcal{T}}\\}$.
For specific classes of correlation functions, [8] presented the closed-form
expressions for the mean and variance of $y$ at new inputs
$\tilde{\boldsymbol{\mathrm{x}}}$ and $\tilde{\boldsymbol{z}}$, conditioned on
$\mathcal{D}$.

We fit univariate GP emulators for the coefficients of the retained basis
vectors obtained for multivariate outputs of component models, and substitute
these parameter estimates in the closed form expression of [8, Theorem 3.1].
Conditioned on the projected ensemble, we obtain a posterior mean and variance
for the output of the composite model for a given input.

##  III  Component models

In this section we provide details for each of the component models of the
decision support framework illustrated in Figure 2 for the case study from a
UK county council considering low carbon technologies to transform its
infrastructure to reach a net-zero carbon target.

$T,E,H$ Heat demand model Energy systems model Gas priceheat demandcost Figure
2: A network of models for energy planning, where $T$ is the surface
temperature, $E$ is the efficiency of the equipment and $H$ is the global
building transmission coefficient.

####  III-1 Heat demand model

Heat demand is calculated based on the degree days statistic (i.e., the sum
over days in which temperature is below a given temperature threshold [15, 16,
17]). We vary the three inputs: surface temperature $(T)$, efficiency of the
equipment $(E)$ and the global building transmission coefficient $(H)$ (which
when multiplied by the degree days gives the heating power required). Surface
temperature takes the form of a time series of annual averages from 2021 to
2050, whereas the other inputs are assumed to be fixed throughout the whole
period. Projections of the annual heating demand are then generated up to
2050. For details on this model we refer to [15].

We derive the domains for the input parameters from the half-hourly energy
data for 2017-2021 provided by the council together with the historical
surface temperature data for 2017-2020.

####  III-2 Energy systems model

This model is based on the widely used OSeMOSYS open source framework [18],
that computes the energy supply mix (in terms of generation capacity and
energy delivery) given the heat demand, fuel prices, and the technologies
connected to the system. The model considers representative days for each of
the four seasons, and daily/seasonal storage technologies which shift demand
from day to night and between seasons, respectively. Based on the half-hourly
demand data provided by the local council, the share of annual heat demand
attributed to each season and day/night is displayed in Table I and specified
in the energy systems model.

TABLE I: Share of annual heat demand (%) | Winter | Spring | Summer | Autumn
---|---|---|---|---
Day | 26.5 | 17.7 | 12.2 | 24.5
Night | 4.66 | 5.11 | 4.12 | 5.14

For this case study, we are interested in output produced by the model of the
total operational cost, a time-series from 2021 to 2050. The model inputs
considered are heating demand and gas price, that take the form of an annual
average time-series from 2021 to 2050 [11].

##  IV Results

###  IV-A Emulation

To construct the decision support system, we first build the GP emulators to
approximate outputs produced by the individual component models in Sections
IV-A1 and IV-A2. In Section IV-B, we outline the process of linking the GP
emulators with mutlivariate output.

####  IV-A1 GP emulator for the heat demand model output

The surface temperature input is represented by one parameter, a shift away
from the central projection on a continuum, where a shift of $\pm 1$
corresponds to the high/low scenario, with intermediate values interpolated
between these [7].

TABLE II: Domain of input parameters for GP emulators for the heat demand model output. Input parameter (unit) | Domain
---|---
Shift parameter in surface temperature | $[-1,1]$
Efficiency of the equipment $(E)$ | $[0.5,1]$
Global building transmission coefficient $(H)$ | $[5,20]$

Table II provides the domains of the input parameters considered for
emulation. We use the maximin distance Latin Hypercube (LHC) to generate a
space-filling design to explore the output behaviour across the input space
and fit the emulator [19]. The output of interest is a vector of annual
heating demands from 2021 to 2050, which is projected onto a low-dimensional
basis by PCA. The first two principal components were found to explain 98% of
the total variance and are retained for our analysis.

We consider coefficients of this basis, denoted by
$c_{i}^{(1)}(\boldsymbol{\mathrm{x}}_{1})$, $i=1,2$, as functions of the input
parameters and construct a GP emulator for each coefficient as in II-A. The
regression function $h(\boldsymbol{\mathrm{x}}_{1})=(1,x_{11},x_{12},x_{13})$
includes a constant, and linear terms in each component of
$\boldsymbol{\mathrm{x}}_{1}$ [1]. The correlation function in the GP is a
squared exponential, and for the fit the RobustGaSP package [20] is used.

To validate the performance of the GP emulators, a test set of size 30 is
used. Figure 3 (top row) presents the cross-validation diagnostics for each
basis coefficient. We plot the emulated values and the model outputs on the
x-axis and y-axis respectively. The black points and error bars represent the
emulator prediction and a two standard deviation prediction interval. The true
model values are green if they lie within two standard deviation prediction
intervals, or red otherwise. We observe that the emulator predictions lie
close to the true values, and the size of the error bars is small for both
basis coefficients, indicating that the emulator is performing well.

Figure 3: Cross-validation plots for emulators for the coefficients of the two
basis vectors for the heat demand model (top) and energy systems model
(bottom).

####  IV-A2 GP emulator for energy systems model output

Similarly, we construct an emulator to approximate the output of the energy
systems model. The heating demand input is replaced by the two basis
coefficients obtained in IV-A1, whereas the gas price input is replaced by a
shift parameter, with high, low and baseline scenarios again represented by
values $\pm 1$ and $0$ of a shift parameter. The output of interest is the
annual operating cost from 2021-50, which is projected on to principal
components, here we use the first two principal components that explain 93% of
the total variance.

The same approach to fitting GP emulators for the retained basis coefficients,
$c_{i}^{(2)}(\cdot),i=1,2$, as a function of two basis coefficients and a
shift parameter for gas price is followed in Figure 3 (bottom row). We observe
that three of the tested model outputs were outside the prediction intervals
for our coefficients of principal components, which is still consistent with
our uncertainty specification.

###  IV-B Linking the GP models

The two emulators for heat demand and the energy systems model are coupled
using the direct calculation method described in Section II-C. For the
graphical model in Figure 2, in Figure 4 we llustrate the relationship between
the first two coefficients of PCA basis vectors obtained from the heat demand
model ensemble, $c_{1}^{(1)}$ and $c_{2}^{(1)}$, and the first two
coefficients of PCA basis vectors obtained from the energy systems model
ensemble, $c_{1}^{(2)}$ and $c_{2}^{(2)}$, given the vectors of global inputs
$\boldsymbol{\mathrm{x}}_{1}$ and $\boldsymbol{\mathrm{x}}_{2}$.

$\boldsymbol{\mathrm{x}}_{1}$$c^{(1)}_{1}$$c^{(1)}_{2}$$c^{(2)}_{1}$$c^{(2)}_{2}$$\boldsymbol{\mathrm{x}}_{2}$
Figure 4: The coefficients for the linked emulators.

Figure 5 presents the cross-validation diagnostics for the linked emulators.
We continue operating with GP emulators for PCA coefficients described in
IV-A1 and IV-A2.

Figure 5: Cross-validation plots for the linked emulators.

Figure 6 shows the projections of annual operational costs. The projected mean
(solid line) and two standard deviation prediction interval (shaded region)
produced by linked emulators (red); and emulation for the energy systems model
only taking full runs of the demand model as inputs (blue). The dashed line is
the simulation run produced by the energy systems model.

The linked emulator allows the consideration of uncertainty from both the heat
demand model and energy systems model; it can be seen in the emulation of only
the energy system model output, the uncertainty of the projection is
underestimated compared to the more comprehensive treatment. A limited
treatment of uncertainty such as the blue series can have significant
consequences in practical decision situations, through consequent projections
being overconfident.

Figure 6:  A simulation run from the energy systems model (dashed line),
projected mean $\pm 2\sigma$ prediction interval produced by the linked
emulators (red) and emulator for the energy systems model only (blue).

##  V Discussion

Simulator models are commonly used to provide decision support for those
managing complex physical systems. We have demonstrated how to analyse two
linked simulator models with multivariate outputs using Gaussian process
emulators, a well-established approach for quantifying uncertainty in computer
models [21]. The presented approach is more efficient compared to performing
Monte Carlo simulation with full probability distributions for the basis
vectors’ coefficients. To our knowledge, this is the first application of such
an approach for decision support within the energy systems domain.

The developed method has been applied to a planning question proposed for a
facility managed by a UK county council, namely a replacement of a gas boiler
with a ground source heat pump. We presented projections for the total
operating costs up to 2050, considering the uncertainties associated with the
environment under which such new systems will need to operate. Under the
carbon zero policy target the cost projections for alternative heating
facilities can be used to support the decision making process in local
councils when various low carbon technologies are considered for long term
planning.

The system model considered for this case study is simple, however the
framework can be applied to larger networks of more complex models within
energy systems planning. Decision support frameworks such as this, which
provide direct expressions for the efficient propagation of uncertainty
between component models are critical for scalability. Further work within
this application domain includes extending the network model to include other
variables, and studying the effect of storage on future projections as well as
smoothing supply-demand fluctuations on operating timescales.

## Acknowledgements

We would like to thank Mark Roberts (Northumberland County Council) for his
constructive and substantive comments.

## References

  * [1] M. Kennedy and A. O’Hagan, “Bayesian calibration of computer models,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 63, no. 3, pp. 425–464, 2001.
  * [2] HM Treasury, “The Aqua Book: guidance on producing quality analysis for government,” _https://www. gov. uk/government/publications/the-aqua-book-guidance-on-producing-qualityanalysis-for-government (Accessed 19-11-2021)_ , 2015.
  * [3] M. Leonelli. and J. Q. Smith, “Bayesian decision support for complex systems with many distributed experts,” _Ann. Oper. Res._ , vol. 235, no. 1, pp. 517–542, 2015.
  * [4] S. Conti, J. Gosling, J. Oakley, and A. O’Hagan, “Gaussian process emulation of dynamic computer codes,” _Biometrika_ , vol. 96, no. 3, pp. 663–676, 2009\.
  * [5] V. Volodina and D. Williamson, “Diagnostics-driven nonstationary emulators using kernel mixtures,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 8, no. 1, pp. 1–26, 2020.
  * [6] D. Williamson and A. Blaker, “Evolving Bayesian emulators for structured chaotic time series, with application to large climate models,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 2, no. 1, pp. 1–28, 2014.
  * [7] A. Wilson, C. Dent, and M. Goldstein, “Quantifying uncertainty in wholesale electricity price projections using bayesian emulation of a generation investment model,” _Sustainable Energy, Grids and Networks_ , vol. 13, pp. 42–55, 2018.
  * [8] D. Ming and S. Guillas, “Linked gaussian process emulation for systems of computer models using matérn kernels and adaptive design,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 9, no. 4, pp. 1615–1642, 2021.
  * [9] F. Sanson, O. Le Maitre, and P. Congedo, “Systems of Gaussian process models for directed chains of solvers,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 352, pp. 32–55, 2019.
  * [10] K. Kyzyurova, J. Berger, and R. Wolpert, “Coupling computer models through linking their statistical emulators,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 6, no. 3, pp. 1151–1171, 2018.
  * [11] “Energy and emissions projections,” https://www.gov.uk/government/collections/energy-and-emissions-projections, 2021, Accessed: 2021-03-29.
  * [12] “Public Sector Decarbonisation Scheme,” https://www.gov.uk/government/publications/public-sector-decarbonisation-scheme-psds, 2021, Accessed: 2021-03-16.
  * [13] R. Haylock and A. O’Hagan, “On inference for outputs of computationally expensive algorithms with uncertainty on the inputs,” _Bayesian statistics_ , vol. 5, pp. 629–637, 1996.
  * [14] D. Higdon, J. Gattiker, B. Williams, and M. Rightley, “Computer model calibration using high-dimensional output,” _Journal of the American Statistical Association_ , vol. 103, no. 482, pp. 570–583, 2008.
  * [15] M. De Rosa, V. Bianco, F. Scarpa, and L. Tagliafico, “Heating and cooling building energy demand evaluation; a simplified model and a modified degree days approach,” _Applied Energy_ , vol. 128, pp. 214–229, 2014.
  * [16] M. Larsen, S. Petrović, A. Radoszynski, R. McKenna, and O. Balyk, “Climate change impacts on trends and extremes in future heating and cooling demands over europe,” _Energy and Buildings_ , vol. 226, p. 110397, 2020.
  * [17] J. Spinoni, J. Vogt, P. Barbosa, A. Dosio, N. McCormick, A. Bigano, and H. Füssel, “Changes of heating and cooling degree-days in Europe from 1981 to 2100,” _International Journal of Climatology_ , vol. 38, pp. e191–e208, 2018.
  * [18] M. Howells, H. Rogner, N. Strachan, C. Heaps, H. Huntington, S. Kypreos, A. Hughes, S. Silveira, J. DeCarolis, and M. Bazillian, “OSeMOSYS: the open source energy modeling system: an introduction to its ethos, structure and development,” _Energy Policy_ , vol. 39, no. 10, pp. 5850–5870, 2011\.
  * [19] M. Morris and T. Mitchell, “Exploratory designs for computational experiments,” _Journal of statistical planning and inference_ , vol. 43, no. 3, pp. 381–402, 1995.
  * [20] M. Gu, J. Palomo, and J. Berger, _RobustGaSP: Robust Gaussian stochastic process emulation_ , 2020, R package version 0.6.1. [Online]. Available: https://CRAN.R-project.org/package=RobustGaSP
  * [21] D. Williamson, M. Goldstein, and A. Blaker, “Fast linked analyses for scenario-based hierarchies,” _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , vol. 61, no. 5, pp. 665–691, 2012.