# A Review of Sequential Decision Making via Simulation

Zhuo Zhang  School of Data Science, the Chinese University of Hong Kong,
Shenzhen, Shenzhen, 518172, China Dan Wang  Department of Industrial
Engineering, School of Mechanical Engineering, Northwestern Polytechnical
University, Xi’an 710072, China Haoxiang Yang  School of Data Science, the
Chinese University of Hong Kong, Shenzhen, Shenzhen, 518172, China Shubin Si
Department of Industrial Engineering, School of Mechanical Engineering,
Northwestern Polytechnical University, Xi’an 710072, China Corresponding
author<EMAIL_ADDRESS>(S. Si).

###### Abstract

Optimization via simulation has been well established to find optimal
solutions and designs in complex systems. However, it still faces modeling and
computational challenges when extended to the multi-stage setting. This survey
reviews the models and methodologies of single-stage optimization via
simulation and multi-stage stochastic programming. These are necessary
theoretical components to push forward the development of sequential decision
making via simulation. We identify the key challenge of sequential decision
making via simulation as the appropriate modeling of the stage-wise value
function, for which we survey the state-of-the-art meta-models and their
potential solution algorithms.

_Keywords_ : optimization under uncertainty; sequential decision process;
multi-stage stochastic program; optimization via simulation.

##  1 Introduction

Optimization via simulation refers to using simulation techniques and
methodologies to optimize complex systems or processes, especially under the
consideration of uncertainty. For many applications, it is difficult to
directly obtain the analytical form of the performance measure function in
terms of decision variables, i.e., the system configuration. The decision
maker makes use of computer-based or physical simulation models to analyze and
evaluate such performance measures under different scenarios and then
identifies the optimal or near-optimal decision variables that result in the
best performance according to predefined criteria, such as expected cost,
equity, reliability, and other utilities. There is usually a realistic limit
for the round of simulation due to a physical budget, for example, the test of
engine design [95], or a computational limit in time and resources [32]. The
objective of the simulation optimization often lies within the following three
categories: 1) to achieve the smallest regret compared to the optimal policy
in the long run, under the assumption that the whole system dynamics are
known; 2) to maximize the probability of choosing the optimal system
configuration under a fixed simulation budget; and 3) to minimize the number
of simulation runs to select the optimal solution with a fixed probabilistic
guarantee.

Optimization via simulation has gained attention in both simulation and
optimization societies because of its broad applicability in many real-world
problems; for example, chemical engineering processes where the fundamental
chemical and physical relationships can be well approximated by simulation
processes [25, 23, 24], transportation and supply chains where simulation
models characterize the flow dynamics and agent behavior [118, 136], drug
design where molecular dynamics simulation assists optimizing the drug
structure and the combination of small compounds [2, 161], manufacturing
processes where complex material flow relationships are captured by simulation
models to help determine optimal machine schedules [130, 155, 103, 102], and
pandemic control where it is natural to use a simulation model in the place of
the differential equations representing pandemic transmission dynamics [117,
120, 160]. A common feature of the applications mentioned above of
optimization via simulation methods is that since it is hard to find a closed-
form representation of many processes in science, engineering, and human
society, researchers resort to characterizing a simulation system to
approximate the reality and hope to make the optimal decision based on the
feedback from the simulation system.

We can formulate a general model of optimization via simulation as follows:

$\displaystyle\min_{x\in{\cal X}}\quad$
$\displaystyle{\mathbb{E}}_{\xi}[f(x,\xi)],$ (1)

where we aim to find the optimal solution $x$ within a feasible region ${\cal
X}$. Even with the same solution $x$, the uncertainty from the environment
$\xi$ will affect the value of the objective function $f$, and such
uncertainty may or may not depend on the choice of $x$. Since we do not know
the realization of the uncertainty $\xi$ when we make the decision, we aim to
optimize the expected value. The expected value operator can be substituted by
other utility functions such as conditional value at risk (CVaR) in risk-
averse optimization [108] and maximization in robust optimization [27].

When the closed form of $f$ is not known and evaluating such a function relies
on simulation outputs, traditional optimization via simulation methods can be
applied. We will detail those methods in Section 2. However, many optimization
via simulation methods do not consider a sequential decision making setup,
i.e., multi-stage, because of the complexity of the convoluted value function
structure without a closed form. Therefore, previous literature reviews mainly
focus on simulation techniques and convergence proofs in a single-stage or a
two-stage setting and do not utilize the potential structure properties of the
value functions within a sequential decision process. We believe it is also
helpful to discuss the multi-stage optimization under uncertainty solution
techniques so that we can better combine them with the simulation methodology
and obtain new algorithms. We start our review on the multi-stage optimization
under uncertainty methods with a multi-stage extension of model (1) as
follows:

$\displaystyle\min_{x_{1}\in{\cal
X}_{1}}\;c_{1}(x_{1})+{\mathbb{E}}_{\xi_{2}}[\min_{x_{2}\in{\cal
X}_{2}(x_{1},\xi_{2})}\;c_{2}(x_{2})+{\mathbb{E}}_{\xi_{3}}[\cdots{\mathbb{E}}_{\xi_{T}}[\min_{x_{T}\in{\cal
X}_{T}(x_{1},\dots,x_{T-1},\xi_{2},\dots,\xi_{T})}c_{T}(x_{T})]]].$ (2)

We can also present model (2) in a recursive way:

$\displaystyle\min_{x_{1}\in{\cal
X}_{1}}\;c_{1}(x_{1})+{\mathbb{E}}_{\xi_{2}}[f_{2}(x_{1},\xi_{2})],$ (3)

where the value function $f_{t}(x_{1},\dots,x_{t-1};\xi_{2},\dots,\xi_{t})$ is
an optimization problem over the future time horizon $\\{t,\dots,T\\}$:

$\displaystyle
f_{t}(x_{1},\dots,x_{t-1};\xi_{2},\dots,\xi_{t})=\min_{x_{t}\in{\cal
X}_{t}(x_{1},\dots,x_{t-1},\xi_{2},\dots,\xi_{t})}c_{t}(x_{t})+{\mathbb{E}}_{\xi_{t+1}}[f_{t+1}(x_{1},\dots,x_{t};\xi_{2},\dots,\xi_{t+1})]$
(4)

At stage $t=1$, the first-stage decision $x_{1}$ is made based on the known
constraint information ${\cal X}_{1}$ and cost function $c_{1}$. However, the
decision $x_{1}$ aims to optimize more than just the first-stage cost but the
expected future cost. At each following stage $t=2,\dots,T$, the uncertainty
$\xi_{t}$ realizes, which helps define the feasible region of the stage-$t$’s
decision and cost function. The multi-stage stochastic program (or its
variants of risk-averse or distributionally robust optimization) often relies
on a known and well-structured closed form of $f$, for example, $f(x,\xi)$ is
a convex linear program in a stochastic programming setting, optimization
methods focus on how to approximate the function $f$. We discuss recent
advances in optimization under uncertainty, especially under the multi-stage
setting, in Section 3.

When we want to extend the sequential decision making process to consider
black-box future value functions $f$, there are some main conflicts between
the current optimization via simulation and sequential decision making. The
sequential decision making process often relies on a rigid framework and
assumptions, which cannot approximate any process to an arbitrary precision.
The optimization via simulation methods focus more on the statistical analysis
and do not fully utilize the structures of the developed optimization
algorithms, e.g., Benders decomposition and dual decomposition algorithms.

Therefore, in this paper, we review the recent progress in optimization via
simulation and sequential decision processes in the hope that we can develop
new methods to combine their advantages and solve complicated multi-stage
problems. We will further discuss potential opportunities to combine the two
research streams via Bayesian optimization and other meta-models in Section 4.

##  2 Review of Traditional Optimization via Simulation Work

Based on the structure and property of the feasible region ${\cal X}$, we may
divide optimization via simulation (OvS) problems into three categories.

  * •

If ${\cal X}$ has a small number of solutions (often less than 100) and the
decision vector $x$ may be numerical or categorical, then we may simulate all
solutions and select the best among them. This problem is known as the
ranking-and-selection (R&S) problem.

  * •

If ${\cal X}$ is a convex subset of $\mathbb{R}^{d}$ and $x$ is a vector of
continuous decision variables, this problem is known as the continuous OvS
(COvS) problem.

  * •

If ${\cal X}$ is a set of $d$-dimensional integers and $x$ is a vector of
discrete-valued variables, this problem is known as the discrete OvS (DOvS)
problem.

###  2.1 Ranking & Selection

There are many reviews available on R&S, such as [15, 66, 89]. Solution
procedures can be categorized as two-stage procedures [42, 133], two-stage
procedures with screening [113], and fully-sequential procedures [88], all of
which can guarantee a desired probability of correct selection. Most
procedures above employ a concept of indifference zone (IZ), where a
difference is considered significant if larger than a pre-specified
indifference-zone parameter. The probability of correct selection guarantee in
the IZ approach corresponds to the probability of selecting the true best,
subject to the condition that the mean of the true best is better than the
mean of all other alternatives by at least the indifference-zone parameter.
Thus, this is based on a worst-case performance metric. This worst-case
approach can provide frequentist guarantees for correct selection but might
require more samples to obtain that guarantee than may be practically
implementable. As such, IZ R&S procedures can be statistically conservative.

There are also R&S procedures based on average-case performance metrics that
sample in a highly sequential manner. The goal is either to maximize evidence
for correct selection subject to a constraint on the sampling budget or to
reach a level of evidence for correct selection with the fewest expected
number of samples. Several typical approaches use an average-case analysis
rather than a worst-case bound of the IZ approach, such as optimal computing
budget allocation (OCBA) proposed by [33] and expected value of information
(EVI) proposed by [34].

As mentioned above, many existing R&S procedures are designed to solve small-
or medium-scale problems. Researchers begin investigating parallel computing
for R&S problems, such as [107, 81, 115, 84, 163, 164]. Meanwhile, some
research problems that expand classical R&S from different perspectives have
emerged, e.g., considering multiple performance measures by multi-objective
R&S [14, 48, 9], taking the input uncertainty into account [53, 157, 46], and
treating the performance measure as a function of the underlying contexts,
which is called R&S with covariates or contextual R&S [148, 52, 149, 37].

###  2.2 Gradient Simulation

For optimization problems with continuous-valued decision variables, the
availability of gradient information can significantly improve the
effectiveness of solution algorithms. However, in the stochastic setting,
since the outputs are themselves random, deriving stochastic gradient
estimators can be challenging.

The stochastic approximation (SA) algorithm is often used to solve
$\mathrm{COvS}$ problems iteratively. Each iteration takes the following form
of update:

$x_{n+1}=\Pi_{{\cal
X}}\left[x_{n}-a_{n}\widehat{\nabla}f\left(x_{n}\right)\right],$

where $x_{n}$ is the solution found at iteration
$n,\widehat{\nabla}f\left(x_{n}\right)$ is an estimate of the gradient $\nabla
f\left(x_{n}\right),\left\\{a_{n}\right\\}$ is a sequence of positive real
numbers called the gain sequence, and $\Pi_{{\cal X}}$ denotes a projection
operator on to the set ${\cal X}$. We can view SA as the stochastic
counterpart to the steepest descent in deterministic optimization. SA is first
introduced in [134] to solve noisy root-finding problems and is later applied
to the setting of stochastic optimization by solving for the zero of the
gradient. Kiefer-Wolfowitz (KW) SA algorithm by [87] addresses the gradient-
free setting.

We can divide the approaches to stochastic gradient estimation into two main
categories, which are indirect and direct.

When the simulation model is treated as a black box, $\nabla f(x)$ may only be
estimated using finite-difference (FD) approximations. Let $e_{i}$ denote the
$i$ th column of a $d\times d$ identity matrix. Then a forward FD estimator of
$\nabla f(x)$ is

$\widehat{\nabla}f(x)=\frac{1}{h}\left[\left(\begin{array}[]{c}Y\left(x+he_{1}\right)\\\
Y\left(x+he_{2}\right)\\\ \vdots\\\
Y\left(x+he_{d}\right)\end{array}\right)-Y(x)\right],$

where $h$ is set as a small positive value. The forward FD estimator is
biased, and the bias is of $O(h)$. When $Y\left(x+he_{i}\right)$ and $Y(x)$
are independent, the variance is of $O\left(h^{-2}\right)$. When the dimension
of $x$ is large, the KW SA algorithm is inefficient. To overcome the
difficulty, Ref. [151] proposes a simultaneous perturbation (SP) FD estimator
called SPSA that uses only two simulation runs to obtain an estimate of
$\nabla f(x)$ regardless of the dimension of $x$.

When the internal structure of the simulation model is known, we may be able
to use this knowledge to design more efficient gradient estimators. Among many
different methods, the most commonly used are perturbation analysis (PA) and
the likelihood ratio/score function (LR/SF) method. Ref. [73] proposes PA
based on the fact that

$\nabla f(x)=\nabla\mathrm{E}[Y(x)]=\mathrm{E}[\nabla Y(x)],$

when the random function $Y(x)$ is stochastically Lipschitz continuous and
differentiable with probability one. Then, if $\nabla Y(x)$ can be observed in
the simulation, PA estimates $\nabla f(x)$ by $\nabla Y(x)$. Several methods
have been proposed for problems where $Y(x)$ is not stochastically Lipschitz
continuous, including the smoothed perturbation analysis in [51] and [67], and
the pathwise method in [76].

Ref. [65] and [132] first introduce the LR/SF method. We can often represent
$Y(x)$ by $Y(x)=h(Z)$, where $Z$ is a vector of random variables generated in
the simulation and $h(Z)$ is the performance measure in which we are
interested. Let $f_{z}(z,x)$ denote the probability density of $Z$. Then,
under some regularity conditions,

$\displaystyle\nabla f(x)$ $\displaystyle=\nabla\mathrm{E}[h(Z)]=\nabla\int
h(z)f_{z}(z,x)dz=\int h(z)\nabla_{x}f_{z}(z,x)dz$ $\displaystyle=\int
h(z)\nabla_{x}\log\left[f_{z}(z,x)\right]f_{z}(z,x)dz=\mathrm{E}\left\\{h(Z)\nabla_{x}\log\left[f_{z}(Z,x)\right]\right\\}.$

If $f_{z}(z,x)$ is known, the LR/SF method estimates $\nabla f(x)$ by
$h(Z)\nabla_{x}\log\left[f_{z}(Z,x)\right]$. Similar to the PA estimator, the
LR/SF estimator is also unbiased and only requires a single simulation run to
compute. Also, it does not require $Y(x)$ to be stochastic Lipschitz
continuous. Therefore, its applicability is wider than the PA estimator.
Despite the benefits of LR/SF, it has poor variance properties, its variance
can be shown to scale linearly with the scale of $Z$ [110].

In recent years, there have been advancements in gradient estimation
techniques. Using FD schemes, it takes $O(d)$ function evaluations to obtain
just one gradient evaluation, which can be expensive for large $d$ even if
parallel computations can be performed. Thus, randomized finite difference
methods (RFD) are proposed to estimate gradients using an arbitrary number of
function value samples and taking an average of directional derivative
estimates along random directions, e.g., [47, 19, 1, 105]. In addition,
several variance reduction techniques for estimators are proposed, such as
control variates [153, 12, 101] and alternative variational objectives [111,
109].

###  2.3 Global Methods

For discrete optimization via simulation, we use a set of general problem-
solving strategies named meta-heuristics. Meta-heuristic algorithms attempt to
locate decent feasible solutions to an optimization problem efficiently. To
this end, they evaluate potential solutions and perform a series of operations
on them to find different, better solutions. We would like to point out that
there are extensive works conducted on global methods in the domains of
machine learning, but we exclusively focus on presenting the commonly employed
global methods in the simulation field due to the scope of this review.

Random search methods involve sampling a set of feasible system designs in
each iteration, conducting simulations of the sampled designs to estimate the
performance of these designs, and then using the simulation results to decide
on what designs should be sampled in the next iteration and on the current
estimate of the optimal system design. Consequently, the class of optimization
algorithms under consideration is broad enough to be applied to solve both
deterministic and stochastic optimization problems with either discrete or
continuous decision parameters (or both). There are many random search
algorithms for DOvS problems, including the stochastic ruler method of [159],
the random search methods of [4], the simulation annealing algorithm of [3],
the stochastic comparison method of [68], the nested partitions algorithms of
[127] and [150], the model reference adaptive search (MRAS) of [79], the
COMPASS algorithms of [77] and [78]. One desirable feature of random search
methods for discrete simulation optimization is that such methods can be shown
to converge almost surely to the set of global optimal solutions of the
underlying (discrete) simulation optimization problem under very general
conditions, see [5, 7]. In addition, random search methods have other
desirable features, including the ability to move rapidly within the feasible
region to identify areas worthy of further investigation, to focus the search
on desirable subsets of the feasible region that have not been extensively
explored, and to adapt the information available in the current optimization
problem; see [6] for a more detailed review for such features.

Another strategy to find good solutions is to use information on the past
search progress and record this information in memory structures.
Metaheuristics that use this strategy are commonly grouped under the umbrella
term tabu search algorithms [62, 63, 64]. A tabu list records the last
encountered solutions and forbids these solutions from being revisited as long
as they are on the list. Alternatively, the tabu list may also record the last
moves that have been made to prevent them from being reversed. The tabu list
(tenure) is constructed specifically for the problem structure, and sometimes,
it might be hard to design such a list to achieve efficiency.

Simulated annealing (SA) is an iterative search method inspired by metal
annealing [90]. The algorithm starts with an initial solution and performs a
stochastic partial search of the state space using appropriate perturbations
and evaluation functions. In each iteration $n$, a random candidate solution
$x^{\prime}_{n}$ is selected from the neighborhood of the current solution
$x_{n-1}$. Uphill moves are accepted with a probability
$\exp{(-\left[f(x^{\prime}_{n})-f(x_{n-1})\right]/T)}$ where $f(\cdot)$ is the
objective function value (to be maximized) of the solution between brackets
and $T$ is an endogenous parameter called temperature. The probability of
acceptance of uphill moves decreases as $T$ decreases. At high temperatures,
the search is almost random, while at low temperatures, the search becomes
almost greedy, i.e., only good moves are accepted [90]. In addition to
requiring that the temperature $T$ be strictly positive for all $n$, most of
the simulated annealing literature assumes that $T_{n}\to 0$ as $n\to\infty$
at a logarithmic rate, e.g. [70].

Genetic algorithms [74], particle swarm optimization [86], and ant colony
algorithms [39] are also widely employed global methods in simulation. These
methods utilize distinct strategies to explore the search space and find
optimal solutions. Genetic algorithms simulate the process of biological
evolution, employing operations such as selection, crossover, and mutation to
explore the search space and find optimal solutions. Particle swarm
optimization algorithms simulate collective behavior observed in bird flocks
or fish schools, continuously adjusting particle positions and velocities to
search for the optimal solution. Ant colony algorithms simulate the foraging
behavior of ants, mimicking ant movement and the release of pheromones in the
search space to search for the optimal solution. The selection of a suitable
method depends on the specific characteristics and requirements of the problem
at hand.

##  3 Review of Traditional Sequential Decision Process Work

The origin of optimization under uncertainty models can be traced back to
Dantzig’s 1955 paper to incorporate random variables as model parameters in
linear programs [36]. Since then, many variants have been developed, such as
stochastic programming [26, 146], robust optimization [17], risk-averse
optimization [139], chance-constrained optimization [31, 119], and
distributionally robust optimization [129]. Many of the papers focus on a two-
stage setting, where the decision-makers optimize the planning decision in the
first stage, observe the uncertainty, and then make the second-stage
adjustment [147, 142, 162, 104]. The models listed above are generally hard to
directly solve via commercial modeling languages and solvers and thus require
transformation, approximation, and decomposition to solve them. It is
impossible to list all the methods to solve optimization under uncertainty
problems, so we list four fundamental ones here: 1) transformation: for
example, robust counterparts derivation for robust optimization models with
row-wise uncertainty [22, 20]; integer programming [106] and nonlinear
programming [121] reformulations for the chance-constrained programs; 2) value
function approximation via inner and outer approximation methods [124, 147,
114]; bounding techniques based on relaxations [80, 55]; 3) decomposition
methods such as L-shaped (Benders) cutting-plane methods [156, 57] and dual
decomposition methods [30, 135]; 4) stochastic approximation that uses
stochastic gradient information for the descent direction [98], where we
essentially consider the stochastic program as a COvS problem but we do have
structural information about the second-stage value function to obtain the
gradient easily.

When we extend the optimization under uncertainty models to a multi-stage
setting, the solution methods become significantly complicated, and some
previously described two-stage transformations no longer apply. For example,
generating a stochastic gradient for a particular stage’s value function is
challenging because the following stages’ value functions have not been
exactly characterized. Only recently Ref. [99] develops a stochastic
approximation algorithm that can generalize for stochastic programs with more
than three stages, but it requires a nested process of sub-gradient
evaluation, which does not scale as the number of stages grows. Therefore, it
is generally hard for the iterative algorithms widely used in nonlinear
optimization and COvS to satisfy the objective of solving multi-stage
problems, and we put our focus on the recent advances in models and solution
methods to characterize the form of value functions. For multi-stage
stochastic programming problems, most solution methods focus on the value
function approximation since the expression
${\mathbb{E}}_{\xi_{t+1}}[f_{t+1}(x_{1},\dots,x_{t};\xi_{2},\dots,\xi_{t+1})]$
does not provide a directly solvable formulation. Ref. [122] first proposes
the stochastic dual dynamic program (SDDP), which iteratively generates a
valid lower approximation at the sampled solutions. For a multi-stage
stochastic program, the expected value functions are approximated from below
by linear cutting-planes and we are able to solve the following problem:

$\displaystyle\underline{f}_{t}(x_{1},\dots,x_{t-1};\xi_{2},\dots,\xi_{t})=\min_{x_{t}\in{\cal
X}_{t}(x_{1},\dots,x_{t-1},\xi_{2},\dots,\xi_{t})}\quad$ $\displaystyle
c_{t}(x_{t})+\theta_{t+1}$ (5a) s.t.
$\displaystyle\theta_{t+1}\geq\pi_{t+1,l}^{\top}x_{t}+v_{t+1,l}\qquad
l=1,2,\dots$ (5b)

The cut coefficients, $\pi_{t+1,\cdot}$ and $v_{t+1,\cdot}$, are obtained from
solving lower level optimization problems $\underline{f}_{t+1}$, i.e.,
_backward pass_ , at the solutions to the sample path problem
$\hat{x}_{1},\dots,\hat{x}_{T}$, simulated and obtained in _forward pass_. If
we assume convexity in the value function of every stage, SDDP guarantees
asymptotic convergence [126, 60], and Ref. [97] provides a thorough discussion
about the complexity of SDDP. Because of the computational viability of SDDP
and its variants, such cutting-plane algorithms have found wide applications
in energy system planning [61]. In recent years, researchers have incorporated
more risk-averse aspects in the multi-stage modeling, including the multi-
stage distributionally robust optimization [125, 44, 10] and the multi-stage
risk-averse optimization [124, 145]. Under a risk-averse setting, the time
consistency issue has been central to guarantee an implementable policy [145,
75, 138], which requires careful modeling to achieve a composite one-step risk
measure structure. For the multi-stage robust optimization, previous works
mainly focus on adaptive decision rules or finite decision space, as the
multi-stage robust optimization model is convex but infinite-dimensional [18,
59]. Ref. [58] develops a new cutting-plane algorithm similar to SDDP but
maintains both upper and lower bounds for the value functions. Another
important direction is the attempt to address the nonconvexity of the value
function. If the value function $f_{t}$ is no longer convex, the linear cut
(5b) may not provide a tight approximation. For the multi-stage stochastic
mixed-integer programs, Ref. [166] proposes Lagrangian cuts, which can be
tight and valid if the state variables $x_{t}$ are binary.

Another important stream of models dealing with multi-stage structured
uncertainty is the Markov decision process. Solution methods consist of exact
dynamic programming techniques, such as value iteration and policy iteration,
and approximate dynamic programming techniques where the value functions are
approximated by a combination of basis functions [85], heuristic functions
[21], or $Q$ functions characterized by machine learning models [29]; see
[128] for a detailed description of dynamic programming and approximate
dynamic programming and [152] for Q-learning and other reinforcement learning
techniques. Recent developments, first proposed in [123], started to integrate
the Markov chain probabilistic structure with a value function iterative
approximation algorithm, such as the policy graph definition of the
uncertainty [40], shown in Figure 1, and the partially observable state in
[41].

Figure 1: An illustration of the Markovian policy graph in [40].

##  4 Towards Sequential Decision Process via Simulation

In many real-world industrial applications, the process is too complicated to
model with high fidelity, and it is usually costly to construct a full-scale
simulation model. In addition, we may need to run such simulation models
multiple times for each system design, i.e., candidate solution, in an
optimization model to evaluate the expected function value. Even with a
successful estimation of the function value (the zeroth order information), we
are not able to directly use it in an optimization algorithm similar to those
listed in Section 3. Therefore, we need approximation models that can
accurately characterize the simulation process, are fast to calculate, and
easily capture the functional relationship between the simulation output and
the design points. When optimizing system design decisions, we can quickly
deploy approximation models in place of the complex simulation model to guide
the design into a close neighborhood of optimality. Such approximation is also
called _meta-model_ or _surrogate model_ , and it plays a significant role in
many engineered systems’ optimization [82, 13, 38].

Many heuristics have been developed for a specific engineering problem, but
there is very scarce literature about the theory considering each stage’s
transition mechanism and objectives as black-box functions. The difficulty
lies in globally approximating the stagewise value functions and embedding the
approximation in the optimization problem such that we can search for the
optimal solution utilizing known optimization algorithms. We list three
different types of meta-models and discuss their advantages and limitations.
Some of the listed meta-models are based on local first- or second-order
approximations, and thus, we need to fit a series of approximations in a
nested optimization scheme, which may not be efficient if the number of stages
is large. Despite this, we believe it is valuable to cover the local meta-
model methodologies, the philosophy of which may inspire the development of
new algorithms.

###  4.1 Response Surface Methodology (RSM)

The _response surface methodology (RSM)_ is a widely applied iterative method
to approximate the input-output simulation process [91, 112]. In each
iteration, RSM constructs a local first-order or second-order polynomial
regression model to approximate the simulation output function locally in a
neighborhood of the incumbent solution. Suppose we denote the value function
to approximate as $f(x)$, e.g., the expected future value function in a multi-
stage stochastic program. Without loss of generality, we can write a second-
order model to approximate $f(x)$ locally:

$g(x)=\beta_{0}+\beta_{1}^{\top}x+x^{\top}B_{2}x,$ (6)

where $\beta_{0},\beta_{1}$ and $B_{2}$ represent coefficients for the
zeroth-, first-, and second-order terms (we could set $B_{2}=0$ to yield a
first-order meta-model). Given a local region to explore, we first select
input design points, $\hat{x}_{1},\hat{x}_{2},\dots$, based on design of
experiment techniques and obtain the estimated expected simulation responses
via multiple rounds of simulations executed at each solution,
$f(\hat{x}_{1}),f(\hat{x}_{2}),\dots$ [94]. The coefficients are fitted at the
simulation design points by minimizing the least square errors between the
estimation $f(\hat{x})$ and the prediction $g(\hat{x})$. The local model
provides a gradient descent, and we can perform one step of the gradient
descent method, with the step size determined by binary search or backtracking
[116, 140]. At the beginning of the RSM procedure, we usually fit a first-
order meta-model to quickly converge to a neighborhood that is likely to
contain the optimum. Towards the end of the RSM process, a second-order
polynomial regression is fitted to capture the functional property more
accurately at the local minimum. We can run multiple parallel searching
processes starting at different input points to avoid being stuck at one local
optimum. Towards the end of the RSM procedure, Ref. [93] discusses the test of
KKT conditions at the obtained design point to check its optimality. RSM can
also be extended to include multiple random responses, selecting one as the
objective and others to formulate constraints [91].

###  4.2 Stochastic Kriging

Stochastic kriging models are another type of general purpose meta-models
[35]; see [92] for a detailed review of kriging models used in simulation and
prediction. It is initially proposed to address the spatial interpolation
problem in the field of mining geology but later used to fit a function with
the sampled data. In the fitting process, the linear or low-order polynomial
regression models tend to have an issue of fitting well locally but generating
a non-robust mapping from the simulation output to the meta-model globally.
The higher-order nonlinear models work well for certain applications, such as
queuing systems, but require domain knowledge and specialized fitting
algorithms and may overfit [8]. The use of kriging in the stochastic
simulation is first introduced in [154]: the key idea of kriging is to assume
that we can calculate all responses with the sum of a deterministic mean
function $f(x)^{\top}\beta$, a random field $M(x)$ that exhibits spatial
correlation, e.g., a Gaussian process, and a noise $\epsilon(x)$:

$Y(x)=f(x)^{\top}\beta+M(x)+\epsilon(x).$ (7)

We observe $n_{j}$ function values $Y_{ji},i=1,\dots,n_{j}$ at a set of
selected design points $\\{\hat{x}_{j}\\}$ for $j=1,\dots,k$. The predictor
$\hat{Y}$ at a new design point $x_{0}$ is a linear function of the mean
simulated values $\bar{Y}_{j}=\sum_{i=1}^{n_{j}}Y_{ji}$:

$\hat{Y}(x_{0})=\lambda_{0}(x_{0})+\sum_{j=1}^{k}\lambda_{j}(x_{0})\bar{Y}_{j}.$
(8)

The coefficients $\lambda_{0},\dots,\lambda_{k}$ correspond to the solution
that minimizes the mean squared error of the predictor assuming the structure
of (7); the fitting process for a Gaussian process $M(x)$ with a power
exponential kernel can be found in [13].

The fitted kriging model can drive the search process for the optimal
solution. Ref. [154] describes how we can fit a kriging model to the
simulation data to search for an optimal design. It involves iteratively
simulating new samples according to the fitted kriging model, re-fitting the
kriging model according to the newly acquired simulated data, and determining
whether KKT conditions are satisfied.

###  4.3 Gaussian Process and Bayesian Optimization

Gaussian processes have been widely used for approximating the nonlinear
functions from dynamical systems because of their readiness to serve as the
prior probability distribution in iterative Bayesian updates; see an
introduction for Gaussian process regression in [16]. A Gaussian process
$M(x)$ can be considered a distribution over functions with a continuous
domain $x\in{\cal X}$, defined by a mean function ${\mathbb{E}}[M(x)]=\mu(x)$
and a covariance function $\Sigma(x_{1},x_{2})=cov(M(x_{1}),M(x_{2}))$ to
characterize the covariance between two random variables obtained from the
Gaussian process at any two points $x_{1}$ and $x_{2}$ in the domain. Gaussian
processes is a global meta-model, which yields a continuously differentiable
function for prediction. The smoothness property of Gaussian processes can
effectively reduce the mean square error compared to other polynomial
regression models in trajectory modeling [143].

Since the conjugate of a normal distribution is still normal, we can update
the Gaussian process estimation by the Bayes law. Suppose we have simulated at
a series of design points $x_{1},\dots,x_{n}$ and obtained their responses
$f(x_{1}),\dots,f(x_{n})$. By proposing a mean function $\mu_{0}(x)$ and a
covariance function (or kernel) $\Sigma_{0}(x,y)$, which should grant a large
positive correlation to two points $x$ and $y$ that are close in the input
space, we yield a prior distribution on $[f(x_{1}),\dots,f(x_{n})]$ as:

$f(x_{1:n})\sim Normal(\mu_{0}(x_{1:n}),\Sigma_{0}(x_{1:n},x_{1:n})),$ (9)

where we follow the notation in [50] to let the subscripts $1:n$ represent the
vectorized values. Suppose in the $n+1$ observation, we would like to acquire
the information at a new design point $x$. We can calculate the conditional
distribution of $f(x)$ based on the observed data:

$\displaystyle f(x)|f(x_{1:n})$ $\displaystyle\sim
Normal(\mu_{n}(x),\sigma_{n}^{2}(x))$ (10a) $\displaystyle\mu_{n}(x)$
$\displaystyle=\Sigma_{0}(x,x_{1:n})){\Sigma_{0}(x_{1:n},x_{1:n}))}^{-1}(f(x_{1:n}-\mu_{0}(x_{1:n}))+\mu_{0}(x)$
(10b) $\displaystyle\sigma_{n}^{2}(x)$
$\displaystyle=\Sigma_{0}(x,x)-\Sigma_{0}(x,x_{1:n})\Sigma_{0}(x_{1:n},x_{1:n})^{-1}\Sigma_{0}(x_{1:n},x).$
(10c)

[16] list a detailed derivation process of this conditional distribution,
i.e., the posterior distribution $f(x)|f(x_{1:n})$ and discussed methods to
select the appropriate mean function and kernel function, e.g., Matérn kernel
and squared exponential kernel. There are many approaches to optimize hyper-
parameters in the kernel function, such as maximum a posteriori (MAP) estimate
[56], log marginal likelihood approach [131], and fully Bayesian approach
[137]. [131] also discuss more efficient algorithms to perform this Bayesian
update step.

Based on the Bayesian update of the Gaussian process, Ref. [96] proposes a
search algorithm, i.e., _Bayesian optimization_ , which is later popularized
by the foundational work of Efficient Global Optimization (EGO) in [83]. Ref.
[28], [144], [50], and many others thoroughly review the Bayesian optimization
approach, stating the situation for which Bayesian optimization is suitable,
including low dimension for the design point $x$, a feasible region with a
simple structure, and continuity. Ref. [50] summarizes the Bayesian
optimization algorithm by the following steps:

  1. Step 1

Initialize with a Gaussian process prior on $f$, including the definition of
the mean and kernel functions $\mu_{0}$ and $\Sigma_{0}$;

  2. Step 2

Observe $f$ at $n_{0}$ design points and set $n=n_{0}$;

  3. Step 3

Update the posterior probability distribution on $f$ based on the observation;

  4. Step 4

Select the maximizer of the acquisition function, $x_{n+1}$, and observe
$f(x_{n+1})$;

  5. Step 5

Increment $n$ and repeat Step 3-4 until the iteration limit.

In each step, the random variable that represents the future value function
$f(x,\xi)$ is approximated by a Gaussian process posterior. The posterior
Gaussian process is then used as the prior input for the next iteration. The
acquisition function guides the search algorithm to the next candidate design
point $x_{n+1}$, where we obtain new samples to update the estimation. The
acquisition function is chosen to balance the exploration and the
exploitation, a concept commonly practiced in multi-armed bandit problems and
reinforcement learning [69]. If maximizing the acquisition function leads to a
new design point that is too close to the best design point among the
evaluated ones, the search algorithm may not be able to escape the local
minimum and explore other potential regions. Common acquisition functions
include the expected improvement [83], knowledge-gradient [49], entropy search
[71], and predictive entropy search [72]. While the expected improvement and
its variants are the most widely used, other alternatives provide validity
under non-standard assumptions.

It requires adaptation of vanilla Bayesian optimization algorithms to run in
parallel and handle noisy measurements and complicated constraints, which are
often necessary for real-world applications and highly related to optimization
under uncertainty. For constraints without noise, [141] define the expected
improvement to be zero if a candidate design point is infeasible. Such
expected constrained improvement is also independently developed by [54]. Ref.
[100] develops a Bayesian expected improvement step for noisy observations and
constraints. The integration step over the posterior of the acquisition
function is achieved via a quasi-Monte Carlo approximation. Ref. [45] proposes
a state-of-the-art scalable Bayesian optimization method that employs trust
regions to confine the search points and transformations for the objective
function and constraints for easier feasibility/optimality detection.

###  4.4 Discussion: Bridging between Simulation Meta-models and Sequential
Decision Making

Most meta-models aim to solve a single-stage optimization problem with a
similar form of model (1), as it is difficult enough to develop a converging
algorithm that simultaneously performs functional estimation and optimum
search while providing a statistical guarantee. Few previous works extend the
meta-models listed in Section 4.1-4.3 to a sequential, or even simply two-
stage, setting. Ref. [11] applies RSM to approximate the objective function of
a two-stage stochastic linear program with recourse, with candidate solutions
selected by central composite design (CCD) and the expected value function
estimated by sample average approximation and Latin hypercube sampling.
However, such a method does not aim to solve the two-stage stochastic program
but uses RSM to perform a sensitivity analysis on the optimal solution.
Recently, Ref. [165] represents the value functions in an approximate dynamic
programming setting via RSM and generalized polynomial chaos models for a
stochastic energy management problem. Ref. [158] develops a two-stage
optimization via simulation method. In the recourse problem, the objective
function is replaced by a local Gaussian process. The expectation of recourse
functions is then approximated by a global Gaussian process. Iterating between
the global-local meta-models, the algorithm can guarantee the convergence to
the optimum for two-stage discrete stochastic programs.

There are still some challenges to extending the current optimization via
simulation methods to a multi-stage setting. Here, we try to list a few points
we believe could require more future research:

  * •

Many of the simulation optimization methods can select multiple design points
and compare them, then finalize the best decision. We disregard the cumulative
objective values at the trial design points as long as the best decision can
be obtained. On the other hand, we only have historical data for data-driven
sequential decision making and need to optimize the current trial. We believe
it is necessary to identify the correct goal before we propose new
optimization methods. Do we have a fixed budget for trial-and-error and only
care about the optimality of the final decision after the sampling process? Or
do we want to find the optimal solution only based on the historical data? Or
do we rather want to optimize the long-run performance where each trial’s
optimal value also counts, but we would also like to reveal more information?
Different applications may find one of those three to be the most suitable.

  * •

The noise in the meta-models’ observations may have a significantly smaller
magnitude than the uncertainty in sequential decision making. This may require
many more samples to obtain an accurate functional evaluation and makes it
hard to achieve any reasonable result when the sample budget or the historical
data is limited.

  * •

It is hard for meta-models to deal with the nonconvexity and feasibility
requirements. The value function is evaluated at every stage by solving an
optimization problem with the future value functions embedded in it. It is
computationally difficult to formulate the future value functions via Gaussian
processes and other meta-models because they most likely will induce
nonconvexity and cannot characterize the feasible region of the current stage
in which we guarantee to find a feasible solution in all future stages for any
uncertainty realization.

  * •

The recursive relationships between stages make the accuracy of meta-model
approximation decay fast. One potential idea to improve the accuracy is to use
machine learning models as surrogates since they can approximate functions
with high precision, especially deep neural networks, but it is still
computationally challenging since they may require a mixed-integer program or
other nonconvex representation, e.g., in Neur2SP [43]. This may also face
issues when the data points are insufficient.

##  5 Conclusions

In this survey, we review the development of methodology in optimization via
simulation and sequential decision processes. As pointed out by many previous
literature, the fast development of computer simulation models makes them more
accurate and applicable for a wide variety of applications, but embedding such
high-fidelity simulation processes in sequential decision making still needs
more effort, especially in handling the nonconvexity and infeasibility. It is
impossible to include every paper with a related topic. Still, we hope this
survey can inspire new research to bridge the gap for sequential decision
making via simulation.

## References

  * [1] A. Akhavan, E. Chzhen, M. Pontil, and A. Tsybakov.  A gradient estimator via l1-randomization for online zero-order optimization with two point feedback.  Advances in Neural Information Processing Systems, 35:7685–7696, 2022.
  * [2] H. Alonso, A. A. Bliznyuk, and J. E. Gready.  Combining docking and molecular dynamic simulations in drug design.  Medicinal Research Reviews, 26(5):531–568, 2006.
  * [3] M. H. Alrefaei and S. Andradóttir.  A simulated annealing algorithm with constant temperature for discrete stochastic optimization.  Management Science, 45(5):748–764, 1999.
  * [4] S. Andradóttir.  A method for discrete stochastic optimization.  Management Science, 41(12):1946–1961, 1995.
  * [5] S. Andradóttir.  A global search method for discrete stochastic optimization.  SIAM Journal on Optimization, 6(2):513–530, 1996.
  * [6] S. Andradóttir.  An overview of simulation optimization via random search.  Handbooks in Operations Research and Management Science, 13:617–631, 2006.
  * [7] S. Andradóttir.  Simulation optimization with countably infinite feasible regions: Efficiency and convergence.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 16(4):357–374, 2006.
  * [8] B. Ankenman, B. L. Nelson, and J. Staum.  Stochastic kriging for simulation metamodeling.  In 2008 Winter Simulation Conference, pages 362–370. IEEE, 2008\.
  * [9] E. A. Applegate, G. Feldman, S. R. Hunter, and R. Pasupathy.  Multi-objective ranking and selection: Optimal sampling laws and tractable approximations via score.  Journal of Simulation, 14(1):21–40, 2020.
  * [10] R. Arora and R. Gao.  Data-driven multistage distributionally robust optimization with nested distance: Time consistency and tractable dynamic reformulations.  Available at Optimization Online, 2022.
  * [11] T. G. Bailey, P. A. Jensen, and D. P. Morton.  Response surface analysis of two-stage stochastic linear programming with recourse.  Naval Research Logistics (NRL), 46(7):753–776, 1999.
  * [12] J. Baker, P. Fearnhead, E. B. Fox, and C. Nemeth.  Control variates for stochastic gradient MCMC.  Statistics and Computing, 29:599–615, 2019.
  * [13] R. R. Barton and M. Meckesheimer.  Metamodel-based simulation optimization.  Handbooks in Operations Research and Management Science, 13:535–574, 2006.
  * [14] D. Batur, L. Wang, and F. F. Choobineh.  Methods for system selection based on sequential mean–variance analysis.  INFORMS Journal on Computing, 30(4):724–738, 2018.
  * [15] R. E. Bechhofer, T. J. Santner, and D. M. Goldsman.  Design and Analysis of Experiments for Statistical Selection, Screening, and Multiple Comparisons.  John Wiley & Sons, New York, 1995.
  * [16] T. Beckers.  An introduction to Gaussian process models.  arXiv preprint arXiv:2102.05497, 2021.
  * [17] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski.  Robust Optimization, volume 28.  Princeton University Press, 2009.
  * [18] A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski.  Adjustable robust solutions of uncertain linear programs.  Mathematical Programming, 99(2):351–376, 2004.
  * [19] A. S. Berahas, R. H. Byrd, and J. Nocedal.  Derivative-free optimization of noisy functions via quasi-Newton methods.  SIAM Journal on Optimization, 29(2):965–993, 2019.
  * [20] D. Bertsimas, D. B. Brown, and C. Caramanis.  Theory and applications of robust optimization.  SIAM Review, 53(3):464–501, 2011.
  * [21] D. Bertsimas and R. Demir.  An approximate dynamic programming approach to multidimensional knapsack problems.  Management Science, 48(4):550–565, 2002.
  * [22] D. Bertsimas and M. Sim.  The price of robustness.  Operations Research, 52(1):35–53, 2004.
  * [23] A. Bhosekar and M. Ierapetritou.  Advances in surrogate based modeling, feasibility analysis, and optimization: A review.  Computers & Chemical Engineering, 108:250–267, 2018.
  * [24] L. T. Biegler.  Recent advances in chemical process optimization.  Chemie Ingenieur Technik, 86(7):943–952, 2014.
  * [25] L. T. Biegler, I. E. Grossmann, and A. W. Westerberg.  Systematic methods for chemical process design.  Prentice Hall, Old Tappan, NJ (United States), 1997.
  * [26] J. R. Birge and F. Louveaux.  Introduction to Stochastic Programming.  Springer Science & Business Media, 2011.
  * [27] I. Bogunovic, J. Scarlett, S. Jegelka, and V. Cevher.  Adversarially robust optimization with gaussian processes.  Advances in Neural Information Processing Systems, 31, 2018.
  * [28] E. Brochu, V. M. Cora, and N. De Freitas.  A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  arXiv preprint arXiv:1012.2599, 2010.
  * [29] L. Busoniu, R. Babuska, B. De Schutter, and D. Ernst.  Reinforcement Learning and Dynamic Programming Using Function Approximators.  CRC press, 2017.
  * [30] C. C. Carøe and R. Schultz.  Dual decomposition in stochastic integer programming.  Operations Research Letters, 24(1-2):37–45, 1999.
  * [31] A. Charnes and W. W. Cooper.  Chance-constrained programming.  Management Science, 6(1):73–79, 1959.
  * [32] C.-H. Chen and L. H. Lee.  Stochastic Simulation Optimization: An Optimal Computing Budget Allocation, volume 1.  World scientific, 2011.
  * [33] C.-H. Chen, J. Lin, E. Yücesan, and S. E. Chick.  Simulation budget allocation for further enhancing the efficiency of ordinal optimization.  Discrete Event Dynamic Systems, 10:251–270, 2000.
  * [34] S. E. Chick and K. Inoue.  New two-stage and sequential procedures for selecting the best simulated system.  Operations Research, 49(5):732–743, 2001.
  * [35] N. Cressie.  The origins of kriging.  Mathematical Geology, 22:239–252, 1990.
  * [36] G. B. Dantzig.  Linear programming under uncertainty.  Management Science, 1(3-4):197–206, 1955.
  * [37] L. Ding, L. J. Hong, H. Shen, and X. Zhang.  Knowledge gradient for selection with covariates: Consistency and computation.  Naval Research Logistics (NRL), 69(3):496–507, 2022.
  * [38] J. V. S. do Amaral, J. A. B. Montevechi, R. de Carvalho Miranda, and W. T. de Sousa Junior.  Metamodel-based simulation optimization: A systematic literature review.  Simulation Modelling Practice and Theory, 114:102403, 2022.
  * [39] M. Dorigo.  Optimization, learning and natural algorithms.  Ph. D. Thesis, Politecnico di Milano, 1992.
  * [40] O. Dowson.  The policy graph decomposition of multistage stochastic programming problems.  Networks, 76(1):3–23, 2020.
  * [41] O. Dowson, D. P. Morton, and B. K. Pagnoncelli.  Partially observable multistage stochastic programming.  Operations Research Letters, 48(4):505–512, 2020.
  * [42] E. J. Dudewicz and S. R. Dalal.  Allocation of observations in ranking and selection with unequal variances.  Sankhyā: The Indian Journal of Statistics, Series B, pages 28–78, 1975.
  * [43] J. Dumouchelle, R. Patel, E. B. Khalil, and M. Bodur.  Neur2sp: Neural two-stage stochastic programming.  arXiv preprint arXiv:2205.12006, 2022.
  * [44] D. Duque and D. P. Morton.  Distributionally robust stochastic dual dynamic programming.  SIAM Journal on Optimization, 30(4):2841–2865, 2020.
  * [45] D. Eriksson and M. Poloczek.  Scalable constrained Bayesian optimization.  In International Conference on Artificial Intelligence and Statistics, pages 730–738. PMLR, 2021.
  * [46] W. Fan, L. J. Hong, and X. Zhang.  Distributionally robust selection of the best.  Management Science, 66(1):190–208, 2020.
  * [47] M. Fazel, R. Ge, S. Kakade, and M. Mesbahi.  Global convergence of policy gradient methods for the linear quadratic regulator.  In International Conference on Machine Learning, pages 1467–1476. PMLR, 2018.
  * [48] G. Feldman and S. R. Hunter.  Score allocations for bi-objective ranking and selection.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 28(1):1–28, 2018.
  * [49] P. Frazier, W. Powell, and S. Dayanik.  The knowledge-gradient policy for correlated normal beliefs.  INFORMS Journal on Computing, 21(4):599–613, 2009.
  * [50] P. I. Frazier.  A tutorial on Bayesian optimization.  arXiv preprint arXiv:1807.02811, 2018.
  * [51] M. Fu, J.-Q. Hu, M. Fu, and J.-Q. Hu.  Conditional Monte Carlo gradient estimation.  Conditional Monte Carlo: Gradient Estimation and Optimization Applications, pages 73–131, 1997.
  * [52] S. Gao, J. Du, and C.-H. Chen.  Selecting the optimal system design under covariates.  In 2019 IEEE 15th International Conference on Automation Science and Engineering (Case), pages 547–552. IEEE, 2019.
  * [53] S. Gao, H. Xiao, E. Zhou, and W. Chen.  Robust ranking and selection with optimal computing budget allocation.  Automatica, 81:30–36, 2017.
  * [54] J. R. Gardner, M. J. Kusner, Z. E. Xu, K. Q. Weinberger, and J. P. Cunningham.  Bayesian optimization with inequality constraints.  In ICML, volume 2014, pages 937–945, 2014.
  * [55] H. Gassmann and W. Ziemba.  A tight upper bound for the expectation of a convex function of a multivariate random variable.  Stochastic Programming 84 Part I, pages 39–53, 1986.
  * [56] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin.  Bayesian Data Analysis.  Chapman and Hall/CRC, 2013.
  * [57] A. M. Geoffrion.  Generalized Benders decomposition.  Journal of optimization theory and applications, 10:237–260, 1972\.
  * [58] A. Georghiou, A. Tsoukalas, and W. Wiesemann.  Robust dual dynamic programming.  Operations Research, 67(3):813–830, 2019.
  * [59] A. Georghiou, W. Wiesemann, and D. Kuhn.  Generalized decision rule approximations for stochastic programming via liftings.  Mathematical Programming, 152:301–338, 2015.
  * [60] P. Girardeau, V. Leclère, and A. B. Philpott.  On the convergence of decomposition methods for multistage stochastic convex programs.  Mathematics of Operations Research, 40(1):130–145, 2014.
  * [61] A. Gjelsvik, B. Mo, and A. Haugstad.  Long-and medium-term operations planning and stochastic modelling in hydro-dominated power systems based on stochastic dual dynamic programming.  Handbook of Power Systems I, pages 33–55, 2010.
  * [62] F. Glover.  Tabu search—part I.  ORSA Journal on Computing, 1(3):190–206, 1989.
  * [63] F. Glover.  Artificial intelligence, heuristic frameworks and tabu search.  Managerial and Decision Economics, 11(5):365–375, 1990.
  * [64] F. Glover and G. A. Kochenberger.  Critical event tabu search for multidimensional knapsack problems.  Meta-heuristics: Theory and Applications, pages 407–427, 1996.
  * [65] P. W. Glynn.  Likelihood ratio gradient estimation for stochastic systems.  Communications of the ACM, 33(10):75–84, 1990.
  * [66] D. Goldsman, B. L. Nelson, and J. Banks.  Comparing systems via simulation.  Handbook of Simulation, pages 273–306, 1998.
  * [67] W.-B. Gong and Y.-C. Ho.  Smoothed (conditional) perturbation analysis of discrete event dynamical systems.  IEEE Transactions on Automatic Control, 32(10):858–866, 1987.
  * [68] W.-B. Gong, Y.-C. Ho, and W. Zhai.  Stochastic comparison algorithm for discrete optimization with estimation.  SIAM Journal on Optimization, 10(2):384–404, 2000.
  * [69] A. K. Gupta, K. G. Smith, and C. E. Shalley.  The interplay between exploration and exploitation.  Academy of management journal, 49(4):693–706, 2006.
  * [70] B. Hajek.  Cooling schedules for optimal annealing.  Mathematics of Operations Research, 13(2):311–329, 1988.
  * [71] P. Hennig and C. J. Schuler.  Entropy search for information-efficient global optimization.  Journal of Machine Learning Research, 13(6), 2012.
  * [72] J. M. Hernández-Lobato, M. W. Hoffman, and Z. Ghahramani.  Predictive entropy search for efficient global optimization of black-box functions.  Advances in Neural Information Processing Systems, 27, 2014.
  * [73] Y. C. Ho and X. Cao.  Perturbation analysis and optimization of queueing networks.  Journal of Optimization Theory and Applications, 40(4):559–582, 1983.
  * [74] J. H. Holland.  Genetic algorithms.  Scientific american, 267(1):66–73, 1992.
  * [75] T. Homem-de Mello and B. K. Pagnoncelli.  Risk aversion in multistage stochastic programming: A modeling and algorithmic perspective.  European Journal of Operational Research, 249(1):188–199, 2016\.
  * [76] L. J. Hong and G. Liu.  Pathwise estimation of probability sensitivities through terminating or steady-state simulations.  Operations Research, 58(2):357–370, 2010.
  * [77] L. J. Hong and B. L. Nelson.  Discrete optimization via simulation using compass.  Operations Research, 54(1):115–129, 2006.
  * [78] L. J. Hong and B. L. Nelson.  A framework for locally convergent random-search algorithms for discrete optimization via simulation.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 17(4):19–es, 2007.
  * [79] J. Hu, M. C. Fu, and S. I. Marcus.  A model reference adaptive search method for global optimization.  Operations Research, 55(3):549–568, 2007.
  * [80] C. Huang, W. T. Ziemba, and A. Ben-Tal.  Bounds on the expectation of a convex function of a random variable: With applications to stochastic programming.  Operations Research, 25(2):315–325, 1977.
  * [81] S. R. Hunter and B. L. Nelson.  Parallel ranking and selection.  In Advances in Modeling and Simulation: Seminal Research from 50 Years of Winter Simulation Conferences, pages 249–275. Springer, 2017.
  * [82] M. S. Hussain, A. A. Javadi, A. Ahangar-Asr, and R. Farmani.  A surrogate model for simulation–optimization of aquifer systems subjected to seawater intrusion.  Journal of Hydrology, 523:542–554, 2015.
  * [83] D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13:455–492, 1998.
  * [84] B. Kamiński and P. Szufel.  On parallel policies for ranking and selection problems.  Journal of Applied Statistics, 45(9):1690–1713, 2018.
  * [85] P. W. Keller, S. Mannor, and D. Precup.  Automatic basis function construction for approximate dynamic programming and reinforcement learning.  In Proceedings of the 23rd international conference on Machine learning, pages 449–456, 2006.
  * [86] J. Kennedy and R. Eberhart.  Particle swarm optimization.  In Proceedings of ICNN’95-international conference on neural networks, volume 4, pages 1942–1948. IEEE, 1995.
  * [87] J. Kiefer and J. Wolfowitz.  Stochastic estimation of the maximum of a regression function.  The Annals of Mathematical Statistics, pages 462–466, 1952.
  * [88] S.-H. Kim and B. L. Nelson.  A fully sequential procedure for indifference-zone selection in simulation.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 11(3):251–273, 2001.
  * [89] S.-H. Kim and B. L. Nelson.  Selecting the best system.  Handbooks in Operations Research and Management Science, 13:501–534, 2006.
  * [90] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi.  Optimization by simulated annealing.  Science, 220(4598):671–680, 1983.
  * [91] J. P. Kleijnen.  Response surface methodology for constrained simulation optimization: An overview.  Simulation Modelling Practice and Theory, 16(1):50–64, 2008.
  * [92] J. P. Kleijnen.  Kriging metamodeling in simulation: A review.  European Journal of Operational Research, 192(3):707–716, 2009\.
  * [93] J. P. Kleijnen.  Response surface methodology.  In Handbook of Simulation Optimization, pages 81–104. Springer, 2014.
  * [94] J. P. Kleijnen.  Design and analysis of simulation experiments: Tutorial.  In Advances in Modeling and Simulation: Seminal Research from 50 Years of Winter Simulation Conferences, pages 135–158. Springer, 2017.
  * [95] H.-M. Koegeler, G. Regner, T. Sams, and K. Gschweitl.  Using simulation and optimization tools to decide engine design concepts.  Technical report, SAE Technical Paper, 2000.
  * [96] H. Kushner.  A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise.  Journal of Basic Engineering, 86(1):97–106, 1964.
  * [97] G. Lan.  Complexity of stochastic dual dynamic programming.  Mathematical Programming, 191(2):717–754, 2022.
  * [98] G. Lan, A. Nemirovski, and A. Shapiro.  Validation analysis of mirror descent stochastic approximation method.  Mathematical Programming, 134(2):425–458, 2012.
  * [99] G. Lan and Z. Zhou.  Dynamic stochastic approximation for multi-stage stochastic optimization.  Mathematical Programming, 187(1-2):487–532, 2021.
  * [100] B. Letham, B. Karrer, G. Ottoni, and E. Bakshy.  Constrained Bayesian optimization with noisy experiments.  arXiv preprint arXiv:1706.07094, 2017.
  * [101] V. Liévin, A. Dittadi, A. Christensen, and O. Winther.  Optimal variance control of the score-function gradient estimator for importance-weighted bounds.  Advances in Neural Information Processing Systems, 33:16591–16602, 2020.
  * [102] J. Liu, C. Li, F. Yang, H. Wan, and R. Uzsoy.  Production planning for semiconductor manufacturing via simulation optimization.  In Proceedings of the 2011 Winter Simulation Conference (WSC), pages 3612–3622. IEEE, 2011.
  * [103] R. Liu, X. Xie, K. Yu, and Q. Hu.  A survey on simulation optimization for the manufacturing system operation.  International Journal of Modelling and Simulation, 38(2):116–127, 2018.
  * [104] X. Liu, S. Küçükyavuz, and J. Luedtke.  Decomposition algorithms for two-stage chance-constrained programs.  Mathematical Programming, 157(1):219–243, 2016.
  * [105] A. Lobanov, B. Alashqar, D. Dvinskikh, and A. Gasnikov.  Gradient-free federated learning methods with $l_{1}$ and $l_{2}$-randomization for non-smooth convex stochastic optimization problems, 2023\.
  * [106] J. Luedtke and S. Ahmed.  A sample approximation approach for optimization with probabilistic constraints.  SIAM Journal on Optimization, 19(2):674–699, 2008.
  * [107] J. Luo, L. J. Hong, B. L. Nelson, and Y. Wu.  Fully sequential procedures for large-scale ranking-and-selection problems in parallel computing environments.  Operations Research, 63(5):1177–1194, 2015.
  * [108] A. Makarova, I. Usmanova, I. Bogunovic, and A. Krause.  Risk-averse heteroscedastic Bayesian optimization.  Advances in Neural Information Processing Systems, 34:17235–17245, 2021.
  * [109] V. Masrani, T. A. Le, and F. Wood.  The thermodynamic variational objective.  Advances in Neural Information Processing Systems, 32, 2019.
  * [110] J. A. McGill, B. A. Ogunnaike, and D. G. Vlachos.  Efficient gradient estimation using finite differencing and likelihood ratios for kinetic Monte Carlo simulations.  Journal of Computational Physics, 231(21):7170–7186, 2012.
  * [111] A. Mnih and D. Rezende.  Variational inference for monte carlo objectives.  In International Conference on Machine Learning, pages 2188–2196. PMLR, 2016.
  * [112] R. H. Myers, D. C. Montgomery, and C. M. Anderson-Cook.  Response surface methodology: process and product optimization using designed experiments.  John Wiley & Sons, 2016.
  * [113] B. L. Nelson, J. Swann, D. Goldsman, and W. Song.  Simple procedures for selecting the best simulated system when the number of alternatives is large.  Operations Research, 49(6):950–963, 2001.
  * [114] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on Optimization, 19(4):1574–1609, 2009.
  * [115] E. C. Ni, D. F. Ciocan, S. G. Henderson, and S. R. Hunter.  Efficient ranking and selection in parallel computing environments.  Operations Research, 65(3):821–836, 2017.
  * [116] J. Nocedal and S. J. Wright.  Numerical Optimization.  Springer, 1999.
  * [117] E. O. Nsoesie, R. J. Beckman, S. Shashaani, K. S. Nagaraj, and M. V. Marathe.  A simulation optimization approach to epidemic forecasting.  PloS One, 8(6):e67164, 2013.
  * [118] J. B. Oliveira, R. S. Lima, and J. A. B. Montevechi.  Perspectives and relationships in supply chain simulation: A systematic literature review.  Simulation Modelling Practice and Theory, 62:166–191, 2016.
  * [119] B. K. Pagnoncelli, S. Ahmed, and A. Shapiro.  Sample average approximation method for chance constrained programming: theory and applications.  Journal of Optimization Theory and Applications, 142(2):399–416, 2009.
  * [120] A. Paleshi, K.-H. Bae, G. Evans, and S. Heragu.  A simulation-based optimization approach for mitigation of pandemic influenza.  IISE Transactions on Healthcare Systems Engineering, 7(2):107–120, 2017.
  * [121] A. Peña-Ordieres, J. R. Luedtke, and A. Wächter.  Solving chance-constrained problems via a smooth sample-based nonlinear approximation.  SIAM Journal on Optimization, 30(3):2221–2250, 2020.
  * [122] M. V. F. Pereira and L. M. V. G. Pinto.  Multi-stage stochastic optimization applied to energy planning.  Mathematical Programming, 52(1-3):359–375, 1991.
  * [123] A. B. Philpott and V. L. de Matos.  Dynamic sampling algorithms for multi-stage stochastic programs with risk aversion.  European Journal of Operational Research, 218(2):470–483, 2012\.
  * [124] A. B. Philpott, V. L. de Matos, and E. Finardi.  On solving multistage stochastic programs with coherent risk measures.  Operations Research, 61(4):957–970, 2013.
  * [125] A. B. Philpott, V. L. de Matos, and L. Kapelevich.  Distributionally robust sddp.  Computational Management Science, 15:431–454, 2018.
  * [126] A. B. Philpott and Z. Guan.  On the convergence of stochastic dual dynamic programming and related methods.  Operations Research Letters, 36(4):450–455, 2008.
  * [127] J. Pichitlamken and B. L. Nelson.  A combined procedure for optimization via simulation.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 13(2):155–179, 2003.
  * [128] W. B. Powell.  Approximate Dynamic Programming: Solving the Curses of Dimensionality, volume 703.  John Wiley & Sons, 2007.
  * [129] H. Rahimian and S. Mehrotra.  Distributionally robust optimization: A review.  arXiv preprint arXiv:1908.05659, 2019.
  * [130] R. V. Ramasesh.  Dynamic job shop scheduling: a survey of simulation research.  Omega, 18(1):43–57, 1990.
  * [131] C. E. Rasmussen, C. K. Williams, et al.  Gaussian Processes for Machine Learning, volume 1.  Springer, 2006.
  * [132] M. I. Reiman and A. Weiss.  Sensitivity analysis for simulations via likelihood ratios.  Operations Research, 37(5):830–844, 1989.
  * [133] Y. Rinott.  On two-stage selection procedures and related probability-inequalities.  Communications in Statistics-Theory and methods, 7(8):799–811, 1978\.
  * [134] H. Robbins and S. Monro.  A stochastic approximation method.  The Annals of Mathematical Statistics, pages 400–407, 1951.
  * [135] R. T. Rockafellar and R. J.-B. Wets.  Scenarios and policy aggregation in optimization under uncertainty.  Mathematics of Operations Research, 16(1):119–147, 1991.
  * [136] J. P. Romero, A. Ibeas, J. L. Moura, J. Benavente, and B. Alonso.  A simulation-optimization approach to design efficient systems of bike-sharing.  Procedia-Social and Behavioral Sciences, 54:646–655, 2012.
  * [137] S. Rossi, M. Heinonen, E. Bonilla, Z. Shen, and M. Filippone.  Sparse Gaussian processes revisited: Bayesian approaches to inducing-variable approximations.  In International Conference on Artificial Intelligence and Statistics, pages 1837–1845. PMLR, 2021.
  * [138] A. Ruszczyński.  Risk-averse dynamic programming for Markov decision processes.  Mathematical Programming, 125:235–261, 2010.
  * [139] A. Ruszczyński.  Advances in risk-averse optimization.  In Theory Driven by Influential Applications, pages 168–190. INFORMS, 2013.
  * [140] M. H. Safizadeh and R. Signorile.  Optimization of simulation via quasi-Newton methods.  ORSA Journal on Computing, 6(4):398–408, 1994.
  * [141] M. Schonlau, W. J. Welch, and D. R. Jones.  Global versus local search in constrained optimization of computer models.  Lecture Notes-Monograph Series, pages 11–25, 1998.
  * [142] R. Schultz, L. Stougie, and M. H. Van Der Vlerk.  Two-stage stochastic integer programming: A survey.  Statistica Neerlandica, 50(3):404–416, 1996.
  * [143] E. Schulz, M. Speekenbrink, and A. Krause.  A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions.  Journal of Mathematical Psychology, 85:1–16, 2018.
  * [144] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas.  Taking the human out of the loop: A review of Bayesian optimization.  Proceedings of the IEEE, 104(1):148–175, 2015.
  * [145] A. Shapiro.  On a time consistency concept in risk averse multistage stochastic programming.  Operations Research Letters, 37(3):143–147, 2009.
  * [146] A. Shapiro, D. Dentcheva, and A. Ruszczyński.  Lectures on Stochastic Programming: Modeling and Theory.  SIAM, 2021.
  * [147] A. Shapiro and T. Homem-de Mello.  A simulation-based approach to two-stage stochastic programming with recourse.  Mathematical Programming, 81(3):301–325, 1998.
  * [148] H. Shen, L. J. Hong, and X. Zhang.  Ranking and selection with covariates.  In 2017 Winter Simulation Conference (WSC), pages 2137–2148. IEEE, 2017.
  * [149] H. Shen, L. J. Hong, and X. Zhang.  Ranking and selection with covariates for personalized decision making.  INFORMS Journal on Computing, 33(4):1500–1519, 2021.
  * [150] L. Shi and S. Olafsson.  Nested partitions method for global optimization.  Operations Research, 48(3):390–407, 2000.
  * [151] J. C. Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  IEEE Transactions on Automatic Control, 37(3):332–341, 1992.
  * [152] R. S. Sutton and A. G. Barto.  Reinforcement Learning: An Introduction.  MIT press, 2018.
  * [153] G. Tucker, A. Mnih, C. J. Maddison, J. Lawson, and J. Sohl-Dickstein.  Rebar: Low-variance, unbiased gradient estimates for discrete latent variable models.  Advances in Neural Information Processing Systems, 30, 2017.
  * [154] W. C. Van Beers and J. P. Kleijnen.  Kriging for interpolation in random simulation.  Journal of the Operational Research Society, 54:255–262, 2003.
  * [155] P. J. Van Laarhoven, E. H. Aarts, and J. K. Lenstra.  Job shop scheduling by simulated annealing.  Operations Research, 40(1):113–125, 1992.
  * [156] R. M. Van Slyke and R. J.-B. Wets.  L-shaped linear programs with applications to optimal control and stochastic programming.  SIAM Journal on Applied Mathematics, 17(4):638–663, 1969.
  * [157] D. Wu and E. Zhou.  Fixed confidence ranking and selection under input uncertainty.  In 2019 Winter Simulation Conference (WSC), pages 3717–3727. IEEE, 2019.
  * [158] W. Xie, Y. Yi, and H. Zheng.  Global-local metamodel assisted two-stage optimization via simulation.  arXiv preprint arXiv:1910.05863, 2019.
  * [159] D. Yan and H. Mukai.  Stochastic discrete optimization.  SIAM Journal on Control and Optimization, 30(3):594–612, 1992.
  * [160] H. Yang, Ö. Sürer, D. Duque, D. P. Morton, B. Singh, S. J. Fox, R. Pasco, K. Pierce, P. Rathouz, V. Valencia, Z. Du, M. Pignone, M. E. Escott, S. I. Adler, S. C. Johnston, and L. A. Meyers.  Design of COVID-19 staged alert systems to ensure healthcare capacity with minimal closures.  Nature Communications, 12(1):3767, 2021.
  * [161] W. Yu and A. D. MacKerell.  Computer-aided drug design methods.  Antibiotics: Methods and Protocols, pages 85–106, 2017.
  * [162] B. Zeng and L. Zhao.  Solving two-stage robust optimization problems using a column-and-constraint generation method.  Operations Research Letters, 41(5):457–461, 2013.
  * [163] Y. Zhong and L. J. Hong.  Knockout-tournament procedures for large-scale ranking and selection in parallel computing environments.  Operations Research, 70(1):432–453, 2022.
  * [164] Y. Zhong, S. Liu, J. Luo, and L. J. Hong.  Speeding up Paulson’s procedure for large-scale problems using parallel computing.  INFORMS Journal on Computing, 34(1):586–606, 2022.
  * [165] Y. Zhuo, J. Zhu, J. Chen, Z. Wang, H. Ye, H. Liu, and M. Liu.  RSM-based approximate dynamic programming for stochastic energy management of power systems.  IEEE Transactions on Power Systems, 2022.
  * [166] J. Zou, S. Ahmed, and X. A. Sun.  Stochastic dual dynamic integer programming.  Mathematical Programming, 175(1-2):461–502, 2019.