###  2.10 Forecasting171717This subsection was written by Fotios Petropoulos.

Forecasting is concerned with the prediction of unknown/future values of one
or multiple variables of interest. If the values of these variables are
collected over time, especially/in particular at regular intervals, the
corresponding problem is referred to as time-series forecasting. The outputs
of forecasting models include point estimates as well as expressions of
uncertainty of such estimates in terms of probabilistic forecasts, prediction
intervals, or path forecasts. Forecasting is applied in a wide range of
applications. In this subsection, we offer an overview of established
forecasting approaches that are useful in social settings (Makridakis et al.,
2020a), such as forecasts produced to support decision making in operations
and supply chain management (§3.12; §3.24, finance (§3.9), energy (§3.19), and
other domains.

Exponential smoothing is one of the most popular families of models for
univariate time-series forecasting. The underlying principle of exponential
smoothing models is that, at every step, the forecast is updated such that the
most recent information is taken into account by exponentially discounting
information from previous periods. The estimates for the exponential smoothing
parameters are based on in-sample fits. The first and simplest exponential
smoothing method, simple (or single) exponential smoothing, was developed by
Brown (1956). This method was able to handle level-only data (no trend nor
seasonal patterns). Soon after, it was extended to handle trended and seasonal
data (Holt, 2004; Winters, 1960). Forty years later, Hyndman et al. (2002)
introduced a fully fledged family of exponential smoothing models that are
represented in a state-space framework. Usually, three states are considered:
level, trend, and seasonality. The way that these three states interact to
produce the final forecast determines the types of trend and seasonality (such
as additive or multiplicative). Exponential smoothing models are fast to
compute and perform well in a wide range of data (Makridakis et al., 2020b),
rendering them ideal benchmarks for forecasting applications. Detailed reviews
of exponential smoothing models are offered by Gardner (2006) and Hyndman et
al. (2008).

Autoregressive integrated moving average (ARIMA) is another very popular
family of univariate forecasting models (for a seminal work on ARIMA, see Box
et al., 1976). In ARIMA, the data are first rendered stationary through
transformations and differencing. The stationary data are then fitted in
linear regression models (see also the next paragraph on regression models) in
which the predictors are either past values of the data (autoregressive terms)
or past errors (moving average terms). ARIMA models are theoretically
appealing as they can depict a wide range of data generation processes. While
manually identifying an optimal ARIMA model can be sometimes challenging,
nowadays automated approaches exist (see, for example, Hyndman and Khandakar,
2008; Franses et al., 2014)

When the variable of interest is known to be affected by other factors (also
called “exogenous variables”), then causal modelling can be applied. In its
simplest form, causal models can be linear or nonlinear regression models that
regress the values of the dependent variable on the values of the independent
variable(s). Apart from the ordinary least squares regression models, other
types of regression models exist, such as the ordinal, logistic, Poisson,
negative binomial regression models as well as the Generalised Linear Models
(GLMs). The dependent variable (variable of interest) is usually continuous,
however specific regression models exist for ordinal or binary dependent
variables, such as the ordinal logistic regression model.” But of course there
are also regression approaches for count data, like Poisson regression or
negative binomial regression (Hilbe’s textbook of the same name is my go-to
reference on this), or more generally Generalized Linear Models (GLMs). I
would assume these to be more relevant to OR than binary or ordinal logistic
regressions.

A common rule for using regression models for forecasting purposes is that the
values of the independent variables are either known or can be predicted, as
is very common in energy forecasting; see Weron (2014) and §3.19.
Transformations of the dependent or independent variables are sometimes
necessary so that assumptions regarding normality of errors and constancy of
the error variance are satisfied (Lago et al., 2021). Another common issue in
regression models is that of multicollinearity between independent variables.
Linear regression models can also be used to produce time-series forecasts
when no exogenous variables are available. In these cases, we can construct
predictors for trend and seasonality and use these predictors as independent
variables to model the time-series patterns. Finally, it is also worth
mentioning that ARIMA models can be extended to ARIMAX models that can include
the effects of exogenous variables, just like autoregressive (AR) models can
be extended to ARX.

Instead of forecasting each time series separately, several approaches exist
in order to forecast time series data as a collection. Multivariate models
(also known as structural models) are designed to model cross-sectional data,
producing forecasts for many variables of interest at the same time. Such
forecasts take into account interactions between all series. A common example
is the vector autoregressive (VAR) models (Sims, 1980; Hasbrouck, 1995).
Another very popular cross-sectional approach is hierarchical forecasting
(Athanasopoulos et al., 2020). Hierarchical forecasting deals with time-series
data that are naturally arranged in hierarchical structures (for example,
product or geographical hierarchies). Forecasts for each node of the hierarchy
are first produced independently using standard univariate forecasting
approaches (such as exponential smoothing or ARIMA); then, forecasts across
the hierarchy are reconciled to achieve coherency (Wickramasuriya et al.,
2019; Hollyman et al., 2021). Hierarchical forecasts offer better accuracy and
are directly relevant for decision makers at multiple levels of an
organisation. A different form of forecasting using multiple series, which is
widely applied in machine-learning methods, is called cross-learning. This
approach implies learning (usually through features; Montero-Manso et al.,
2020; Wang et al., 2022c) from other series to be able to predict the variable
of interest. Compared to other cross-sectional approaches, cross-learning
requires access to a set of “reference” data which, though, do not have to be
concurrent to the target data.

Given the plethora of available modelling options, we need ways to help us
decide on the best approach for the target data. Two popular approaches for
model selection are information criteria and cross-validation. Information
criteria select the best model amongst a pool of candidate models based on how
well the in-sample forecasts fit the actual data (model fit), penalising at
the same time for model complexity (Occam’s razor). Information criteria are
fast to compute and widely applied, mostly due to their implementations in
open-source forecasting packages (Hyndman and Khandakar, 2008). Cross-
validation is based on the comparison of the out-of-sample performance between
different models. To achieve this, the available data are split into
“training” and “validation” data. The validation follows a rolling-origin
process, where the forecasts of the candidate models are compared for multiple
forecast origins (Tashman, 2000; Bergmeir and Benítez, 2012). A more recent
approach to forecast selection is based on the concept of representativeness
(Petropoulos and Siemsen, 2022). Out-of-sample forecasts with higher
representativeness to the past data patterns are preferred to ones with lower
representativeness. Regardless of how one selects between forecasts and
models, the values of the selection criteria can also be used to combine
forecasts (Kolassa, 2011). In fact, multiple studies have shown that combining
forecasts, using equal or unequal weights, can significantly boost the
forecasting performance of individual models (Bates and Granger, 1969;
Nowotarski et al., 2016; Wang et al., 2022d). Claeskens et al. (2016) offer a
possible explanation on why the performance of forecast combinations is better
than that of the individual forecasts.

Apart from statistical, algorithmic and computational approaches, the
forecasting process can also be infused by judgement (see, also, §2.2 and
§2.20). It is not unusual for forecasts to be directly produced in a
judgemental way, without the support of any systematic approaches. Research
suggests that such forecasts suffer from several biases (Lawrence et al.,
2006). However, managers may sometimes have a unique appreciation of the
situation, one that the hard data cannot communicate through models. In such
cases, systematic approaches to elicit expert knowledge include the Delphi
method (Rowe and Wright, 1999), structured analogies (Green and Armstrong,
2007), prediction markets (Wolfers and Zitzewitz, 2004), and interaction
groups (Van de Ven and Delbeco, 1971); see also Graefe and Armstrong (2011)
and Nikolopoulos et al. (2015) for a comparison between these approaches.
Apart from producing forecasts directly, judgement may also be used to adjust
the formal/statistical forecasts. Judgemental interventions and their efficacy
have been well-studied in the literature (see, for example: Fildes et al.,
2009; Petropoulos et al., 2016; Fildes et al., 2019). The main takeaways are:
(i) negative adjustments are generally more beneficial than positive ones;
(ii) larger adjustments should be preferred to smaller ones; and (iii) the use
of feedback and support will limit and improve the role of such judgemental
adjustments. Finally, managerial judgement may be applied in other stages of
the forecasting process, such as judgementally selecting between statistical
models (Petropoulos et al., 2018; De Baets and Harvey, 2020) or setting their
(hyper)parameters.

Forecasts produced in previous periods need to be evaluated once the
corresponding actual values become available. Through feedback, forecast
evaluation allows analysts to improve the forecasting process and, thus,
forecasting performance. The main rule of forecast evaluation is that
performance should be measured on data that were not used to fit the models or
produce the forecasts. Solely measuring the in-sample performance will
inevitably lead to over-fitting and the use of complex forecasting models.
There exist a wide array of evaluation metrics. Some of them are suitable for
measuring the accuracy or bias of the point forecasts, while others focus on
how well the uncertainty around the forecasts is estimated. In the former
category, popular metrics are the mean absolute error (MAE), the root mean
squared error (RMSE), the mean absolute percentage error (MAPE, which is very
popular in practice) and the mean absolute scaled error (MASE, which is
theoretically more elegant and popular in academia). It should be noted that
Kolassa (2020) showed that different error metrics are minimised by different
(point) forecasts and that it makes little sense to evaluate one point
forecast using multiple KPIs. For detailed overviews of forecasting metrics
for point forecasts and their proper use, the reader is referred to Hyndman
and Koehler (2006), Davydenko and Fildes (2013), and Koutsandreas et al.
(2022). In the latter category, a popular metric is the interval score (IS),
which is a proper scoring rule and considers both the calibration and the
sharpness of the prediction intervals, as well as the pinball score, the
continuous ranked probability score (CRPS), and the energy score. Gneiting and
Raftery (2007) offer a review of (strictly) proper scoring rules. Finally, we
should mention that nowadays it is common to go beyond strict forecasting
performance and measure the performance of the forecasts on the utility (Hong
et al., 2020; Yardley and Petropoulos, 2021).

For a detailed encyclopedic overview of the forecasting field, both in terms
of theory and practice, we refer the reader to the work of Petropoulos et al.
(2022); a live version of this encyclopedia is available at forecasting-
encyclopedia.com. Hyndman and Athanasopoulos (2021) and Ord et al. (2017) have
written comprehensive textbooks on forecasting and its applications. Notable
open-source packages with implementations of the most popular forecasting
models include the forecast (Hyndman et al., 2022) and smooth (Svetunkov,
2022) packages for R statistical software.

###  2.11 Game theory181818This subsection was written by Georges Zaccour.

Game theory is a branch of mathematics that studies strategic interactions
between decision makers, called players. Strategic interactions means that a
player’s payoff depends not only on her own decision (action or choice), but
also on the decisions made by the other players. The book by von Neumann and
Morgenstern (1944) is often considered the starting date of game theory,
though some of its roots can be traced back to much earlier. Games can be
classified along a series of features. In a static game, each player acts only
once, whereas in a dynamic game, interactions are repeated over time. In a
one-person game, the decision maker plays against a nonstrategic (or dummy)
player, often referred to as “nature”, whose action is the outcome of a
probabilistic event with a fixed (known) distribution. Two-player games focus
on one-on-one interactions. Duopolistic competition and management-union
negotiations are situations that can be modelled as two-person games.
Extending the model to $n>2$ players is conceptually easy but may be
computationally challenging because each player needs to determine all the
possible sequences of actions and reactions for all players. When the number
of interacting players is very large, e.g., an economy with many small agents,
the analysis shifts from individual-level decisions to understanding the
group’s behavioural dynamics. An illustration of this is traffic congestion:
when an agent attempts to minimise her travel time on a route from A to B, her
travel speed depends on the traffic density on that route. What matters is the
number of drivers, not their identity. Population and evolutionary games
(Hofbauer and Sigmund, 1998; Cressman, 2003; Sandholm, 2010) and mean-field
games (Huang et al., 2003, 2006, 2007; Lasry and Lions, 2006a, b, 2007; Gomes
and Saúde, 2014) are branches of game theory that study situations with large
numbers of players.

A game can be defined in three forms, namely, in strategic, extensive, or in
coalitional form. To formulate a one-shot game in strategic form, we have to
specify (i) the set of players and, for each player, (ii) the set of actions,
and (iii) a payoff function measuring the desirability of the game’s possible
outcomes, which depends on the actions chosen by all players. The set of
actions can be finite, e.g., to bid on a contract or not, or continuous, e.g.,
the amount bid. If the players intervene more than once in the game, then we
should additionally define (iv) the order of play, (v) the information
acquired by the players over time (stages), and (vi) whether or not nature is
involved in the game.

In a one-shot game, an action (move) and a strategy mean the same thing. In
games where players intervene more than once, the two concepts no longer
coincide. A strategy is then a decision rule that associates a player’s action
with the information available to her at the time she selects her move. So an
action, e.g., spending advertising dollars, is a result of the strategy. The
word strategy comes from Greek (strategia) and has a military sense. An army
general’s main task is to design a plan that takes into account (adapts to)
all possible contingencies. This is precisely the meaning of strategy in game
theory. Whether in war, business or politics, it is never wise to allow
yourself to be surprised by the enemy. This does not imply that a winning
strategy always exists. Sometimes we must be content with a draw or even a
reasonable loss.

One-shot games are a useful representation of strategic interactions when the
past and the future are irrelevant to the analysis. However, if today’s
decisions also affect future outcomes and are dependent on past moves, then a
dynamic game is needed. In a repeated game, the agents play the same game in
each round, that is, the set of actions and the payoff structures are the same
in all stages (Mertens et al., 2015). The number of stages can be finite or
infinite, and this distinction matters in terms of achievable outcomes. In a
stochastic game, the transition between states depends on the players’ actions
(Shapley, 1953; Mertens and Neyman, 1981; Jaśkiewicz and Nowak, 2018a, b). In
a multistage game, the players share the control of a discrete-time dynamic
system (state equations) observed over stages (Başar and Olsder, 1999;
Engwerda, 2005; Krawczyk and Petkov, 2018). Their choice of control levels,
e.g., investments in production capacity or advertising, affects the evolution
of the state variables (e.g., production capacity, reputation of the firm), as
well as current payoffs. Differential games are continuous-time counterparts
of multistage games (Isaacs, 1975; Başar et al., 2018).

Information plays an important role in any decision process. In a game, the
information structure refers to what the players know about the game and its
history when they choose an action. A player has complete information if she
knows who the players are, which set of actions is available to each one, what
each player’s information structure is, and what the players’ possible
outcomes can be. Otherwise, the player has incomplete information. If, for
instance, competing firms do not know their rivals’ production costs, then the
game is an incomplete-information game. The game can also have perfect or
imperfect information. Roughly speaking, in a game of perfect information,
each player knows the other players’ moves when she chooses her own action, as
in, e.g., chess or a manufacturer-retailer game where the upstream player
first announces a product’s wholesale price, and then the downstream player
reacts by selecting the retail price. The archetype of an imperfect-
information game is the prisoner’s dilemma, where (in the original story) the
players have to simultaneously choose between confessing or denying a crime. A
Cournot oligopoly, where each firm chooses its own production level without
knowing its competitors’ choices, is another instance of an imperfect-
information game.

The outcome of a game depends on the players’ behaviour. In a noncooperative
game, e.g., R&D competition to develop a vaccine, each player optimises her
own payoff, whereas in a cooperative game, the players seek a collectively
optimal solution. For instance, the members of a supply chain could agree to
coordinate their strategies to maximise the chain’s total profit. The
fundamental solution concept in a noncooperative game is the Nash equilibrium
(Nash, 1950b, 1951). Let $I=\left\\{1,\ldots,n\right\\}$ be the set of
players, $S_{i}$ the set of strategies of player $i\in I$, and let her payoff
function be given by $g_{i}\left(s_{1},\ldots,s_{n}\right)$
$:\prod\limits_{i\in I}S_{i}\rightarrow\mathbb{R}$, where
$s=\left(s_{1},\ldots,s_{n}\right)$. Assuming the players are maximisers, the
strategy profile $s^{N}=\left(s_{1}^{N},\ldots,s_{n}^{N}\right)$ is a Nash
equilibrium if $g_{i}\left(s_{1}^{N},\ldots,s_{n}^{N}\right)\geq$
$g_{i}\left(s_{1}^{N},\ldots
s_{i-1}^{N},s_{i},s_{i+1}^{N}\ldots,s_{n}^{N}\right)$ for all $s_{i}\in S_{i}$
and all $i\in I$. At an equilibrium, no player has an interest in deviating
unilaterally to any other admissible strategy. Put differently, if all other
players stick to their equilibrium values, then player $i$ does not regret
implementing her equilibrium value too, which is obtained by best-replying to
the choice of the others. That is, $s_{i}^{N}=\arg\max_{s_{i}\in
S_{i}}g_{i}\left(s_{1}^{N},\ldots
s_{i-1}^{N},s_{i},s_{i+1}^{N}\ldots,s_{n}^{N}\right).$ A Nash equilibrium does
not always exist, and there may be multiple equilibria, raising the question
of which one to select (Selten, 1975). Existence and uniqueness conditions for
Nash equilibrium typically rely on fixed-point theorems. If the game is one of
incomplete information, then the solution concept is a Bayesian Nash
equilibrium (Harsanyi, 1967, 1968a, 1968b). Another noncooperative equilibrium
solution concept, which predates the Nash equilibrium, is the Stackelberg
equilibrium, introduced in a two-player framework by von Stackelberg (1934).
There is a hierarchy in decision-making between the two players: the leader
first announces her action, and next the follower makes a decision that takes
the leader’s action as given. Before announcing her action, the leader would
of course anticipate the follower’s response and selects the action that gives
her the most favourable outcome. The framework has been extended to several
followers and leaders (Sherali, 1984).

In a cooperative game, the players coordinate their strategies in view of
optimising a collective outcome, e.g., a weighted sum of their payoffs, and
must agree on how to share the dividend of their cooperation (Moulin, 1988;
Owen, 1995). Different solution concepts have been proposed, each based on
some desirable properties, typically stated as axioms, such as fairness,
uniqueness of allocation, and stability of cooperation. The most-used
solutions in applications are the core (Gillies, 1953), and the Shapley value
(Shapley, 1953). In any solution, the set of acceptable allocations only
includes those that are individually rational. Individual rationality means
that a player will agree to cooperate only if she can get a better outcome in
the cooperative agreement than she would by acting alone. In a dynamic
cooperative game, the agreement must specify, at the outset, the decisions
that must be implemented by each player throughout the planning horizon. One
concern in such games is the durability of the agreement over time. Clearly,
it is rational for a player to leave the agreement at an intermediate time if
she can achieve a better outcome. The literature on dynamic games has followed
two streams in its quest to sustain cooperation over time, namely, building
cooperative equilibria or defining time-consistent solutions. Through the
implementation of some (punishing) strategies, the first stream seeks to make
the cooperative solution an equilibrium of an associated noncooperative game.
If this is achieved, then the result will be at once collectively optimal and
stable, as no player will find it optimal to deviate unilaterally from the
equilibrium. See Osborne and Rubinstein (1994) for repeated games, Dutta
(1995) and Parilina and Zaccour (2015) for different types of stochastic
games; and Haurie and Tolwinski (1985), Tolwinski et al. (1986), and Haurie
and Pohjola (1987) for multistage and differential games. The second stream
looks for time-consistent solutions, which are achieved by allocating the
cooperative payoffs over time in such a way that, along the cooperative state
trajectory, no player will find it optimal to switch to her noncooperative
strategies. The idea was initiated in Petrosjan (1977) and has since been
further developed (see Yeung and Petrosyan, 2018; Petrosyan and Zaccour,
2018).

Game theory has found applications in biology, economics, engineering,
management, Operational Research, and political and social sciences.

###  2.12 Graphs and networks191919This subsection was written by Ivana
Ljubić.

Graphs and networks are used to represent interactions, connections or
relationships between objects. In network optimisation problems, numerical
attributes representing features such as costs, weights or capacities are
assigned to objects (also called _vertices_) or to connections between them.
If connections are directed, we refer to them as _arcs_ , otherwise we call
them _edges_. Given an input graph with $n$ vertices and $m$ arcs (or edges),
the goal is to find a subgraph that exhibits desired properties (described by
a given set of constraints) and that optimises the given objective function
(usually measured as the sum of edge or vertex “weights” of the solution’s
subgraph). In the following, we focus on some of the most fundamental and most
studied problems in network optimisation.

The _shortest path problem_ in arc-weighted graphs, for example, seeks to find
a least costly path from the given source vertex $s$ to the given target $t$.
When the arc costs are non-negative, one can use the algorithm of Dijkstra
(1959), the efficient implementation of which uses Fibonacci heaps and runs in
$\mathcal{O}(m+n\log n)$ time. For graphs with possible negative arc costs, in
$\mathcal{O}(mn)$ time the Bellman-Ford algorithm either finds the shortest
path from $s$ to all other vertices, or it proves that such a path does not
exist due to the presence of a negative cost cycle reachable from $s$. The
shortest path algorithms are explained in many textbooks, see e.g., Cormen et
al. (2022); Kleinberg and Tardos (2006); Schrijver (2003); Williamson (2019).

In the _maximum flow problem_ (MF), in a given network with arc capacities, we
want to send as much flow as possible from the given source $s$ to the given
sink $t$ without violating the arc capacities. The problem was motivated by
the conflict between East and West during the Cold War (Schrijver, 2002). Ford
and Fulkerson (1957) develop the first exact algorithm that searches for
augmenting paths in the _residual network_. Their fundamental result, known as
the _max-flow/min-cut theorem_ states that the maximum flow passing from the
source to the sink is equal to the total capacity of the arcs in a _minimum
cut_ , i.e., the network that indicates how much more flow is allowed in each
arc., which is a subset of arcs of the smallest total capacity, the removal of
which disconnects the source from the sink. The same result using the duality
theory of LPs is given in Dantzig and Fulkerson (1955). The famous results
from graph theory such as Menger’s theorem, König-Egeváry theorem, or Hall’s
theorem, follow from the max-flow/min-cut theorem (Ford and Fulkerson, 1962).
The method of Ford and Fulkerson (1957) is pseudo-polynomial when arc
capacities are integral, however it may fail to find the optimal solution and
need not terminate if some of the arc capacities are irrational (Ford and
Fulkerson, 1962). An algorithm that overcomes this issue was independently
discovered in the 1970s by Edmonds and Karp (1972) and Dinic (1970), see also
Dinitz (2006). Augmenting the flow along shortest paths (that is, along the
paths with fewest edges) guarantees a polynomial-time complexity. Instead of
augmenting the flow along a single augmenting path as in Edmonds and Karp
(1972), the algorithm of Dinic (1970) finds all shortest augmenting paths in a
single phase. Another stream of MF algorithms exploits the _preflow_ idea of
Karzanov (1974) in which the vertices are overloaded with the excess flow
(i.e., more incoming flow than the outgoing flow is allowed). Subsequent
improvements are obtained in the following years. An important breakthrough is
achieved by Goldberg & Tarjan with the introduction of _push-relabel_
algorithms (Goldberg and Tarjan, 1988). A _pseudoflow_ algorithm for the
maximum flow is introduced by Hochbaum (2008) and it is later improved in
Hochbaum and Orlin (2013). The recent implementation by Goldberg et al. (2015)
is competitive with the Boykov and Kolmogorov (2004) method and the pseudoflow
approach. Further historical details and a more complete list of references
can be found in Ahuja et al. (1993); Dinitz (2006); Goldberg and Tarjan
(2014); Williamson (2019). Currently, the best strongly polynomial bounds are
obtained by Orlin (2013) and King et al. (1994). However, new and improved MF
algorithms continue to be discovered. The most recent trends use the idea of
electrical flows for obtaining faster (exact or approximate) algorithms, see,
e.g., Chapter 8 of Williamson (2019).

In the _minimum cost flow problem_ (MCF), for each arc of the graph, a cost is
incurred per unit of flow that traverses it. The goal is to send units of a
good that reside at one or more _supply_ vertices to some other _demand_
vertices, without violating the given arc capacities at minimum possible cost.
Edmonds and Karp (1972) introduce the _scaling technique_ for the MCF. The
technique is later improved by Orlin (1993). The algorithms of Vygen (2002)
and Orlin (1993) have the best-known strongly polynomial complexity bounds for
the MCF. Kovács (2015) provides a comprehensive literature overview and gives
an experimental evaluation of MCF algorithms based on _network-simplex,
scaling_ or _cycle-cancelling_ techniques. The MCF is treated in detail in
many textbooks, Ahuja et al. (1993); Korte and Vygen (2008); Williamson
(2019). One of the important results is the _integrality of flow_ property: if
all demands/supplies and arc capacities are integers, then there exists an
optimal MCF solution with integer flow on each arc. The result follows from
the _totally unimodular_ property of the constraint matrix when the MCF is
modelled as a linear program.

In the _minimum cut problem_ (MC), one searches for a proper subset of
vertices $S$ of a given arc-capacitated graph, such that the total capacity of
arcs leaving $S$ is minimised. For directed graphs, the algorithm of Hao and
Orlin (1994) is based on MF calculations between chosen pairs/subsets of
vertices and exploits the push-relabel ideas. For undirected graphs, the
_Gomory–Hu tree_ , which is a weighted tree that represents the minimum
$s$-$t$ cuts for every $s$-$t$ pair in the graph, is introduced in Gomory and
Hu (1961). This tree is constructed after $n-1$ MF computations, and a simpler
procedure has been later given by Gusfield (1990). The algorithm of Padberg
and Rinaldi (1990a) improves the ideas of Gomory and Hu (1961) and is widely
used within branch-and-cut schemes for solving the travelling salesperson
problem (TSP) and related problems. The _maximum adjacency ordering_ together
with Fibonacci heaps is used in Nagamochi et al. (1994). Randomised approaches
can be found in Karger and Stein (1996); Karger (2000). The method of Karger
(2000) is de-randomised by Li (2021). Practical performance of some of these
algorithms is evaluated in Chekuri et al. (1997); Jünger et al. (2000). For
additional and more recent references, see the book by Williamson (2019).

The problems mentioned so far all belong to the class $\mathcal{P}$, however
most of the network optimisation problems that are relevant for practical
applications are $\mathcal{NP}$-hard. We highlight two of them that serve as
drivers for discovering new algorithms and methodologies that can be easily
adapted to other difficult optimisation problems.

Given an undirected graph with non-negative edge costs, the _Steiner tree
problem in graphs_ (STP) asks for finding a subtree that interconnects a given
set of vertices (referred to as _terminals_) at minimum cost. Two special
cases can be solved in polynomial time: when all vertices are terminals (the
minimum spanning tree problem), or when there are only two terminals (the
shortest path problem). In general, however, the decision version of the STP
is $\mathcal{NP}$-complete (Karp, 1972b). Older surveys covering developments
of first MIP formulations, Lagrangian relaxations, branch-and-bound methods
and heuristics can be found in Maculan (1987); Winter (1987). The research on
the STP was marked by polyhedral studies in the 1990s (Goemans, 1994; Chopra
and Rao, 1994). Exact solution methods for the STP are based on a
sophisticated combination of: _reduction techniques_ (Gamrath et al., 2017;
Rehfeldt and Koch, 2021), _dual and primal heuristics_ (Pajor et al., 2018)
embedded within branch-and-cut or branch-and-bound frameworks, see (Polzin,
2003; Vahdati Daneshmand, 2003; Polzin and Vahdati Daneshmand, 2009; Gamrath
et al., 2017; Fischetti et al., 2017a). Currently best approximation ratio for
the STP is 1.39 (Goemans et al., 2012). A comprehensive survey of the results
obtained in the last three decades is given by Ljubić (2021). State-of-the-art
computational techniques for the STP are due to Rehfeldt (2021).

The _Travelling salesperson problem_ (TSP) aims at finding the answer to the
following question: If a travelling salesperson wishes to visit all $n$ cities
from a given list exactly once, and then return to the home city, what is the
cheapest route they need to take? For the history of the problem, see
Applegate et al. (2011) and the book by Cook (2011). Since 1954, when Dantzig
et al. (1954) found a provably optimal solution for a 49-city problem
instance, many important improvements in the development of exact methods have
been achieved202020www.math.uwaterloo.ca/tsp/history/milestone.html. Facet-
defining inequalities are investigated in Padberg and Rinaldi (1990b); Jünger
et al. (1995). MIP formulations, including the famous _subtour-elimination
constraints_ model by Dantzig et al. (1954), are compared in Padberg and Sung
(1991). Branch-and-cut methods are developed in Applegate et al. (2011);
Jünger et al. (1995); Padberg and Rinaldi (1991). For the most recent overview
on approximation algorithms for the TSP see Traub (2020). Helsgaun
(2000)212121akira.ruc.dk/~keld/research/LKH/ provides an efficient
implementation of the _$k$ -opt heuristic_ of Lin and Kernighan (1973). Cook
et al. (2021) extend the algorithm of Helsgaun (2000) to deal with additional
constraints in routing applications and win the Amazon Last Mile Routing
Challenge in 2021. The TSP solver
Concorde222222www.math.uwaterloo.ca/tsp/concorde (Applegate et al., 2011)
incorporates best algorithmic ideas from the past 60 years of research on the
topic. By combining techniques of Helsgaun (2000) and Applegate et al. (2011),
instances with millions of vertices can be solved to within 1% of optimality,
see e.g., TSP solutions on graphs with up to 1.33 billion of
vertices232323www.math.uwaterloo.ca/tsp/star.

###  2.13 Heuristics242424This subsection was written by Ceyda Oğuz and İstenç
Tarhan.

Etymologically meaning to find/discover, heuristics make use of previous
experience and intuition to solve a problem. A heuristic algorithm is designed
to solve a problem in a shorter time than exact methods, by using different
techniques ranging from simple greedy rules to complex structures, which could
be dependent on the problem characteristics; however it does not guarantee to
find the optimal solution (§2.4; §2.9). Heuristics have been used in the
operational research area extensively with respect to the applications (see,
for example, §3.12, §3.14, §3.15, and §3.32). In this subsection, we review
the methods employed in the development of heuristics.

Classifications and strategies provided in the literature guide us for the
methods employed in heuristics. Below we provide a thorough classification and
explain briefly the basic methods used under each class.

Induction, being the simplest method to be applied with an analogy to the
mathematical induction, is to solve the original complex problem by extending
the results and insights obtained from small and simpler versions of the
problem (Silver et al., 1980; Silver, 2004; Laguna and Martí, 2013).

Restriction methods primarily focus on explicitly eliminating some parts of
the solution space so that the problem will be solved given a restricted set
of solutions (Silver et al., 1980; Zanakis et al., 1989; Silver, 2004; Laguna
and Martí, 2013). One way of doing this is to identify common attributes of
the optimal solution and search among the solutions having these attributes
only (Glover, 1977). Another restriction can be applied by eliminating
infeasible solutions considering a combination of decision variables which
dictates incompatible values. Beam Search (Morton and Pentico, 1993) is a good
example of this class of heuristics which works with a truncated tree
structure using strategies similar to a branch-and-bound algorithm (§2.4). The
trimming of the tree is utilised by a parameter called beam width to indicate
how many nodes to have at every level of the tree.

Heuristics using decomposition/partitioning method employ different approaches
to divide the problem into smaller and tractable parts, solve these parts
separately and combine their solutions to give the solution to the original
problem (Foulds, 1983; Zanakis et al., 1989; Silver, 2004; Laguna and Martí,
2013). The methods used to divide and then combine the solutions are usually
dictated by the nature of the problem. For example, Hierarchical Planning
proposed by Hax and Meal (1973) considers the organisational level breakdown
and the output of one decomposed problem becomes the input for the other.
Rolling horizon also falls under this category (Stadtler, 2003). A problem
with a sequence of decisions that span a long planning horizon is solved by
dividing the planning horizon into smaller planning intervals. The problem
with these small planning intervals is solved continually by fixing the
decisions for the first time period and moving into the next time period to
solve the next problem. Another approach takes the characteristics of the
input data into account and divides the problem such that each part includes
only tractable amount of data. For example, data showing clusters of
geographically close customers is suitable for this type of partitioning. The
decomposed problems are solved independently, and their solutions are combined
with a certain rule. Divide and Conquer algorithm heuristically clusters
vertices on a given graph, generates a smaller graph for each cluster and
solves the original problem for each cluster independently (Akhmedov et al.,
2016). Decomposition can be made based on an element of the problem, for
example solving a logistics problem after dividing it into parts per vehicle.
Other decomposition approaches benefit from the structure of the mathematical
model developed for the problem. Examples of this sort are Lagrangian
Relaxation (Fisher, 1981), in which complicating constraints are lifted to the
objective function with a penalty, and Benders Decomposition (Benders, 1962;
Rahmaniani et al., 2017), in which once complicating variables are fixed, the
remaining problem can be divided into problems to be solved independently.

Approximation methods focus on the mathematical models and utilise different
strategies to make the problem tractable which results in a reduced size of
the problem (Silver et al., 1980; Silver, 2004). One strategy widely used is
the aggregation over variables or stages. Another common strategy is to modify
the variables, the objective function or the constraints of the mathematical
model in different ways, such as converting discrete variables into continuous
variables, using a linear objective function instead of the non-linear
objective function,, linearising nonlinear constraints, and either eliminating
or weakening some of the constraints (Glover, 1977). Kernel Search (Angelelli
et al., 2010), which combines relaxation with decomposition over the decision
variables, demonstrates that a heuristic may use more than one class of
methods in its design.

Constructive heuristics start from an empty solution and build a complete
solution by adding an element of the problem following a rule at every step,
such as the nearest neighbour algorithm (Bellmore and Nemhauser, 1968) for the
travelling salesman problem. Usually constructive heuristics are of greedy
nature by making the decision for local optimum in every step. These
algorithms can be enhanced by adding a look-ahead mechanism that is by
estimating the future effects of a decision rather than just the current
effect to avoid pitfalls of being greedy.

Improvement heuristics start with a complete solution and improve it by
modifying one or more elements of the solution in every iteration until a
predetermined stopping condition is achieved. Improvement heuristics in their
simplest form utilise a local search which is defined over a neighbourhood
structure to express how the moves are performed from one iteration to the
next. $k$-opt is an example of this sort which replaces $k$ elements of a
solution with another set of $k$ elements in every step if it is beneficial
(Lin and Kernighan, 1973). The parameter $k$ determines the size of the local
search and implicitly applies the restriction method discussed above. A
neighbourhood is defined by a set of solutions which are reachable form the
current solution. A local search is performed by moving from the current
solution to another solution in the neighbourhood of it (next solution).
Selection of the next solution is done by accepting either the one among
random choices that improves the objective function value first (random
descent if it is a minimisation problem) or the one resulting in the best
objective function, i.e., the local optimum, with respect to that
neighbourhood (steepest descent for a minimisation problem). This simple
structure focuses on the local information (exploitation of the accumulated
search experience) and is known as intensification (Glover, 1990). While it
will be useful if the structure of the problem is appropriate, it may result
in not good enough solutions otherwise. Hence, the improvement heuristic will
benefit if it can explore other parts of the solution space, which is known as
diversification (Glover, 1990). Two immediate strategies to be employed are
either to start the search from different initial solutions and choose among
the final solutions obtained (multi-start algorithms) or to allow moving to
worse solutions if this direction will provide a better path for the future
selections (hill-climbing strategy for a minimisation problem).

Even though metaheuristics (Glover, 1986) are improvement methods, since they
advance notably, considering them as a separate class is worthwhile.
Metaheuristics utilise a local search together with intensification and
diversification mechanisms and aim at eliminating the problem-dependent and
domain-specific nature of other heuristics. Simulated Annealing (Kirkpatrick
et al., 1983) is one of the most popular metaheuristics which uses a single
solution in its local search with a random descent and utilises hill-climbing
strategy for diversification. Tabu Search (Glover, 1986) is an example of
deterministic metaheuristic working with a single solution throughout the
search. It explicitly uses history of search in both intensification and
diversification mechanisms. Genetic Algorithm (Holland, 1975) is another
popular metaheuristic comprising of random components for intensification and
diversification but working with a set of solutions during the search.
Variable Neighbourhood Search (Hansen and Mladenović, 1999) is an excellent
example of a design in which diversification is provided by systematically
changing neighbourhood structures.

Matheuristics are heuristic approaches that exploit exact approaches (and
their complementary strengths) without guaranteeing to find the optimal
solutions. While matheuristics are designed with different strategies, we
summarise the most widely used three strategies.

Those matheuristics which are originally exact approaches yet are implemented
heuristically are overlapping with what is described under restriction and
decomposition/partitioning methods in this subsection. Apart from those
overlapping works, in the context of dynamic programming, the corridor method
constructs neighbourhoods as corridors around the state trajectory of the
incumbent solution (Sniedovich and Viß, 2006). Defined (preferably large)
neighbourhoods can be searched with exact approaches. Dynasearch algorithm
uses dynamic programming to search an exponential size neighbourhood stemming
from compound moves in polynomial time (Congram et al., 2002).

Another group of matheuristics benefits from multiple exact models
collectively within a heuristic mechanism. Tarhan and Oğuz (2022) decompose
the scheduling planning horizon into a set of buckets, solve a time-indexed
model to generate a restricted model for each bucket and solve the restricted
models sequentially to construct a complete feasible solution. Della Croce et
al. (2014) solve a restricted time-indexed model and a model with positional
variables iteratively to search the neighbourhood of the incumbent solution.
Solyalı and Süral (2022) propose a matheuristic algorithm by sequentially
solving different mixed integer linear programs.

Third strategy is to incorporate exact models into different components of the
heuristics. This approach may have several variations. First version includes
those matheuristics having a constant interaction between heuristics and
mathematical programming models. Manerba and Mansini (2014) use the Variable
Neighbourhood Search to decide which variables to fix in their fix-and-
optimise algorithm. Archetti et al. (2015) use different integer programming
models in both the intensification and the diversification phases of their
Tabu Search algorithm to improve the objective function value and/or restore
feasibility. Adouani et al. (2022) apply exact and heuristic approaches
respectively to change the value of so-called upper and lower level variables
in the neighbourhood search. Other variations include matheuristics that
sequentially call heuristics and the models; e.g., exact approaches following
heuristics for post-optimisation (Pillac et al., 2013), exact approaches
generating the initial solutions from heuristics (Macrina et al., 2019), exact
approaches supporting heuristics at their both beginning and end to provide an
initial solution and to improve the final solution, respectively (Archetti et
al., 2017).

We refer the reader for a detailed overview of heuristics to the works of
Müller-Merbach (1981) and Silver (2004), of metaheuristics to the work of Blum
and Roli (2003), and of matheuristics to the work of Boschetti and Maniezzo
(2022). The most recent book by Martí et al. (2018) on heuristics is another
invaluable resource. Finally, the progress of metaheuristics is discussed by
Swan et al. (2022). This work provides a critical analysis of the current
state of metaheuristics by focusing on cultural and technical barriers.

For the future studies in the area of heuristics, new techniques and powerful
mechanisms could be derived from practical problems to address complex systems
of today’s world. Another contribution can be to explore and integrate
applications of artificial intelligence to deal with large scale data.
Matheuristics are especially often applied for single-objective problems and
accordingly, their implementation for multi-objective optimisation is a
promising future research direction. For practical purposes, such as to be
used within commercial solvers, it is also worthwhile to develop generic
matheuristic frameworks that can address specific classes of optimisation
problems. Parallel computing (i.e., parallel solution of mathematical models)
and integration with machine learning (to, for example, manage the interaction
with mathematical models and heuristics) are some other invaluable research
directions for matheuristics.

###  2.14 Linear programming252525This subsection was written by Jean-Marie
Bourjolly.

Linear programming (LP) offers a framework for modelling the problem of
extremising a linear economic function under a set of linear inequality
constraints. Solving such models can be approached algebraically as well as
geometrically: finding an extreme point of a polyhedron at which a given
economic function is maximised or minimised. Since its inception in 1947 by
Dantzig, the simplex method has been the standard algorithm for solving linear
programs. A precursor, unbeknownst then to Dantzig, was a set of ideas exposed
by Fourier in 1826 and 1827, and partly rediscovered by Motzkin in 1936, hence
the now famous Fourier-Motzkin elimination method (Dantzig, 1963, p. 84–85;
Schrijver, 1998, p. 155–157) that solves a set of linear inequalities by
sequentially eliminating variables, at the cost though of exponentially
increasing the number of constraints.

But since the 1930s, several researchers had been making a headway. Working
independently from one another, they had grappled with specific problems:
balancing the distribution of revenue (output) with the distribution of
outlays (input) in the economic activity of a whole country (Leontief, 1936);
general economic equilibrium (von Neumann, 1945); production planning
(Kantorovich, 1960); transportation planning (Hitchcock, 1941; Kantorovitch,
1958; Koopmans, 1949); deployment planning and logistics (Dantzig, 1991).
Dantzig (1982) said he had been “fascinated” by Leontief’s interindustry
input-output model and wanted to generalise it by considering many alternative
activities. He also credits von Neumann with the duality theory of linear
programming, which parallels the work the latter did with Morgenstern on the
theory of games.

A linear program can always be expressed (in standard form) as {minimise $cx$,
subject to $Ax=b,x\geq 0$}, where $x$ is an $n$-vector of decision variables,
$A$ is an $m$ by $n$ constraint matrix that somehow weighs the variables, $b$
is an $m$-vector that puts limits on the possible values of $x$, and $cx$ is
an economic function, called the objective function, that measures the quality
of a given solution $x$. It is customary to assume, without loss of
generality, that matrix $A$ is of rank $m$ and that $m$ is smaller than $n$
(see, e.g., Papadimitriou and Steiglitz, 1982). Since
$\text{min}(cx)=-\text{max}(-cx)$, one can minimise a “cost” as well maximise
a “profit”. Linear programs come in pairs: {minimise $cx$, subject to
$Ax=b,x\geq 0$} and {maximise $yb$, subject to $yA\leq c$}. The former is
called the primal, the latter is the dual. The duality theorem of linear
programming has been proved to be equivalent to Farkas’s lemma that was
published in 1902 (see, e.g., Dantzig, 1963). This implies that finding an
optimal solution to a linear program is equivalent to finding a feasible
solution to a system made of the primal and dual constraints with the
additional inequality $cx\leq yb$.

With the introduction of duality, we now have three algorithms for solving
LPs: the (primal) simplex method that maintains (primal) feasibility
throughout and tries to achieve optimality; the dual simplex method that
maintains dual feasibility and moves toward primal feasibility; and the
primal-dual algorithm that starts with a feasible solution to the dual and
keeps improving it by solving an associated restricted primal. The primal-dual
algorithm is the favoured simplex tool for solving most network flow problems,
for instance, the famous algorithm of Ford and Fulkerson (1962) for maximum
flow. The dual simplex method together with column generation may come in
handy when the number of constraints is huge in comparison with the number of
variables (Desaulniers et al., 2006).

A set of linear inequalities defines a convex polyhedron $P$. Therefore, since
the objective function is linear, there are only three possibilities: no
feasible solution (if and only if $P$ is empty); exactly one optimal solution,
located at some extreme point of $P$; infinitely many optimal solutions,
located at the points of a face of $P$ of dimension 1 or more, including its
extreme points. The simplex method moves sequentially along the edges of $P$
from one extreme point to another. Algebraically, it moves from one set of $m$
linearly independent columns, called a basis, to another. Each basis induces a
basic solution defined by setting to zero all the $n-m$ variables that do not
correspond to its columns. A basic solution is feasible if all its components
are non-negative. The move from one basis to another goes as follows: one
column is dropped and replaced by a new one. This exchange, called a pivot,
follows a set of rules for choosing the column that enters the basis and the
one that exits. It is such that, barring degeneracy, the objective function
decreases strictly in value at each pivot.

Degeneracy is rooted in the fact that an extreme point of $P$ may correspond
to several bases. The algebraic expression of this defectiveness is a basic
feasible solution with more than $n-m$ zero components. This occurs when the
number of hyperplanes intersecting at an extreme point is greater than the
minimum necessary to define it. (Think of the tip of a pyramid that has a
square base.) Pivoting in the presence of degeneracy may cause the simplex
method to cycle. Several schemes have been devised to avoid cycling by
carefully choosing the entering and leaving columns. Bland’s rule, considered
as both simple and elegant, has been widely adopted (Bland, 1977). As for
finding an initial solution, if the problem is feasible, this can be done by
introducing artificial non-negative variables that one then tries to drive
down to zero.

Evidence shows that the simplex method is very fast in practice (see Shamir,
1987), but Klee and Minty (1972) designed an LP for which it must visit each
one of the $2^{n}$ or so extreme points, which proved that it is not “good” in
the sense of Edmonds (1965b). A “good algorithm” having been defined as one
for which the worst-case complexity is polynomial with respect to the
dimension of any instance, an important open question became “Is LP in
$\mathcal{P}$?”. Khachiyan answered by the affirmative in 1979 when he adapted
to the specific case of linear programming a known approach in convex
optimisation that had been contributed to by several Soviet mathematicians
(see Gács and Lovász, 1981; Bland et al., 1981; Chvátal, 1983). The argument
goes as follows: given an LP, start with an ellipsoid that is big enough to
contain the set $S$ of feasible solutions if it is not empty. At each
iteration, check whether the centre of the ellipsoid is a solution. If it is
not, there is a hyperplane $H$ separating it from $S$. Cut the ellipsoid in
half by the hyperplane parallel to $H$ that goes through the centre. Then
determine the smallest ellipsoid that contains the half-ellipsoid where one is
trying to locate $S$, and repeat. Stop either with a solution (located at a
centre) or with an ellipsoid that is too small to contain $S$. This is an
important theoretical result (see Grötschel et al., 1981), but with very
little practical use as far as solving actual LPs goes.

The same cannot be said, however, of the interior point algorithm introduced
by Karmarkar (1984), in which the moves happen strictly inside the set of
feasible solutions instead of taking place on the envelope. Indeed,
Karmarkar’s algorithm is polynomial and often competitive with the simplex
method. It assumes a canonical form for linear programming in which the
variables are constrained to $Ax=0,x\geq 0$ and $x\in
S={x:x_{1}+x_{2}+\dots+x_{n}=1}$; it further assumes, without loss of
generality, that the point $e/n=(1/n,1/n,\dots,1/n)$ is feasible and that the
minimum value of the objective function $cx$ is zero. As it seeks to stay away
from the envelope of the solution polyhedron, the algorithm builds a sequence
of strictly feasible solutions, i.e., that have strictly positive components,
and makes a repetitive use of $e/n$. The gist of the algorithm is the
following: given a strictly feasible solution $x^{k}$, one can define a simple
bijective scaling function $f$ that maps $S$ onto itself so that $x^{k}$ is
mapped onto $e/n$, away from the envelope, and so that $f$ has the following
property: if, for any variable $x$, $f(x)$ is strictly feasible in the “new”
space, then so too is $x$ in the initial space. In the “new” space, the
gradient of the transformed objective function is projected on the null space
of the transformed matrix $A$ augmented of a row of 1s, to account for $S$. If
$p$ denotes that projection, one then moves in the direction of $-p$, i.e., in
the direction of the steepest descent, while feasibility is maintained. The
algorithm stops at a point $y^{k+1}$ before reaching the envelope of the
feasible region. That point is transformed to $x^{k+1}$ by $f^{-1}$, and this
is repeated with a new scaling bijection (see Strang, 1987; Goldfarb and Todd,
1989; Fang and Puthenpura, 1993; Winston and Venkataramanan, 2003). Important
links between Karmarkar’s algorithm and the ellipsoid method have been pointed
out (Todd, 1988; Ye, 1987).

Integer Linear Programming (ILP), i.e., linear programs in which the variables
are restricted to being integer-valued, is arguably the most challenging and
beautiful expression of LP. Unfortunately, whereas LP is in $\cal{P}$, ILP is
not, unless $\cal{P}=\cal{NP}$ (Karp, 1972a). However, there are classes of
LPs for which, if there exists a solution at all, an integer solution is
guaranteed without having to make it a requirement. This is the case, e.g., of
most network flow models. And there are classes of ILPs for which any extreme
point of the polyhedron of integer solutions can be obtained by “shaving off”
non-integer extreme points of the outer polyhedron of real-valued solutions
with hyperplanes the number of which is bounded by a polynomial in the
dimension of the instance (Edmonds, 1965a, b; Grötschel et al., 1981, 1988;
Cook et al., 1998). Furthermore, tackling $\cal{NP}$-complete problems has
benefited greatly from this approach (see, e.g., Dantzig et al., 1954;
Applegate et al., 2007).

###  2.15 Integer programming262626This subsection was written by Adam N.
Letchford and Andrea Lodi.

_Mixed-integer programming_ (MIP) is an $\cal{NP}$-hard generalisation of
linear programming (LP; §2.14), in which some or all of the variables are
required to take whole-number values. Way back in the late 1950s, it was
already realised that a wide variety of important practical problems could be
modelled as MIPs (Dantzig, 1960; Markowitz and Manne, 1957). Of course, at the
time, there were no good algorithms, or indeed computers, to enable one to
solve MIPs from real-life applications. Since then, however, dramatic progress
has been made in theory, algorithms and software. Indeed, it is now possible
to solve many real-life MIPs to proven optimality (or at least near-
optimality) on a laptop. In this subsection, we review the main developments
in this area. For more details, we refer the reader to the textbooks by Chen
et al. (2011) and Conforti et al. (2014).

In 1958, Gomory (1958) developed the first finitely-convergent exact algorithm
for pure IPs (i.e., MIPs in which all variables are restricted to whole-number
values). His method was based on _cutting planes_ , i.e. additional linear
constraints which cut off fractional LP solutions. Shortly after, Land and
Doig (1960) invented the _branch-and-bound_ method, in which a sequence of LP
relaxations is embedded within a tree structure. A few years later, Balas
(1965) devised a simpler branch-and-bound algorithm, for pure 0-1 LPs, which
did not rely on LPs at all.

In the 1960s and 1970s, researchers invested considerable effort into deriving
“deep” cutting planes. This led to the discovery of Gomory mixed-integer cuts
(Gomory, 1960), corner polyhedra (Gomory, 1969), intersection cuts (Balas,
1971), Chvátal-Gomory cuts (Chvátal, 1973), disjunctive cuts (Balas, 1979;
Owen, 1973), and cuts derived from a study of the so-called knapsack polytope
(Balas, 1975; Wolsey, 1975). These topics are still being studied to this day
(see, e.g., Conforti et al., 2014; Cornuéjols, 2008).

In 1980, Balas and Martin (1980) developed a general-purpose heuristic for 0-1
LPs, called “pivot-and-complement”. This initiated a line of work on so-called
“primal heuristics”, which also continues to this day. We will mention this
again below.

A major step forward occurred in 1983, with the publication of an award-
winning paper by Crowder et al. (1983). Basically, they did the following
before running branch-and-bound: (i) “pre-process” the formulation in order to
make the LP relaxation stronger, (ii) automatically generate knapsack cuts to
further improve the relaxation, (iii) run a simple primal heuristic in order
to obtain a feasible integer solution early on, and (iv) permanently fix some
variables to 0 or 1 based on reduced-cost arguments. In this way, they were
able to solve ten real-life 0-1 LPs that had previously been regarded as
unsolvable. The largest of these instances had 2756 variables and 756
constraints, a phenomenal achievement at the time. The approach of Crowder et
al. (1983) is now called “cut-and-branch”.

Around the same time, there were several major theoretical advances, such as
the proof of the “polynomial equivalence of separation and optimisation”
(Grötschel et al., 1981) and the development of a polynomial-time algorithm
for pure IPs with a fixed number of variables (Lenstra, 1983). For details, we
recommend Schrijver (1986).

Coming back to a more practical perspective, several improvements were made to
the basic cut-and-branch scheme in the 1980s and 1990s. For brevity, we just
mention some highlights. Several authors proposed more powerful pre-processing
procedures (e.g., Dietrich et al., 1993; Hoffman and Padberg, 1991;
Savelsbergh, 1994). Gu et al. (1998) developed more effective algorithms for
generating knapsack cuts. Researchers also began to study cutting planes for
_mixed_ 0-1 LPs (e.g., Padberg et al., 1984; Van Roy and Wolsey, 1986), which
eventually led to effective cut-and-branch algorithms for such problems (e.g.,
Van Roy and Wolsey, 1987).

The next milestone was the invention of _branch-and-cut_ by Padberg and
Rinaldi (1987). In branch-and-cut, one has the option of generating cutting
planes at _any_ node of the branch-and-bound tree, rather than only at the
root node (as in cut-and-branch). Although this is a fairly simple idea,
Padberg and Rinaldi added several ingredients to turn it into a highly
effective tool. For example, (i) care is taken to ensure that cutting planes
generated at one node of the tree remain valid at all other nodes, (ii)
whenever a cutting plane is generated, it is stored in a so-called “cut pool”,
(iii) when visiting a new node of the tree, one can check the cut pool to see
if it contains any useful cuts, (iv) one uses a heuristic rule to decide when
to stop cutting and start branching at any given node.

Several developments in the 1990s are also worth mentioning. First, there were
some interesting works on methods to construct “hierarchies” of relaxations
for 0-1 and mixed 0-1 LPs (e.g., Balas et al., 1993; Lovász and Schrijver,
1991; Sherali and Adams, 1990). The method in Balas et al. (1993), called
_lift-and-project_ , turned out to be useful when embedded within a branch-
and-cut algorithm for mixed 0-1 LPs (Balas et al., 1996a). Shortly after that,
Balas et al. (1996b) obtained good results using Gomory mixed-integer cuts
instead. This last result was a big surprise: up to then, researchers had
thought that Gomory cuts were of theoretical interest only.

By the end of the 1990s, researchers were routinely solving real-life MIPs
with thousands of variables and hundreds of constraints to proven optimality.
Of course, MIP in general is ${\cal NP}$-hard, so one could not expect to
solve _all_ instances so quickly. Indeed, Cornuéjols and Dawande (1999) found
a family of 0-1 LPs, called “market split” problems, which proved to be
especially challenging for branch-and-cut. This led to the development of a
new class of specific algorithms called _basis reduction_ methods, see, e.g.,
Aardal et al. (2000).

In the period 2000-2010, there was a flurry of impressive works concerned with
primal heuristics for MIP. For brevity, we mention just a few examples.
Fischetti and Lodi (2003) devised a method called _local branching_ , which is
essentially a form of neighbourhood search in which the neighbourhoods – being
of exponential size – are searched by solving auxiliary MIPs. Shortly after,
Danna et al. (2005) presented _relaxation-induced neighbourhood search_ or
RINS, which solves a series of small MIPs to search for integer solutions that
are “close” to the solution of the LP relaxation. Both local branching and
RINS are improving heuristics, i.e., the neighbourhoods are defined with
respect to a reference feasible solution to be improved. Remarkably, they
solve auxiliary MIPs by simply calling a MIP solver in a black-box fashion
(with work limits), thus witnessing the maturity of the field. In the same
year of RINS, Fischetti et al. (2005) introduced the _feasibility pump_ ,
which is highly effective for MIPs where even finding a feasible solution is
challenging.

The development of the branch-and-cut technology has been so impressive that
many of the above-mentioned developments have been incorporated in software
packages. This includes major commercial packages, such as CPLEX, Gurobi and
FICO Xpress, and non-commercial ones that are free to academics, such as SCIP.
We remark that this continual development in algorithms and software has been
greatly enhanced by the creation and constant maintenance of MIPLIB, a library
of MIP instances on which all new methods are now routinely tested (see Bixby
et al., 1992; Gleixner et al., 2021).

We end this section by briefly mentioning three other areas of constant
development. First, there has been great progress on _decomposition_
approaches to MIPs that have special structure, with _branch-and-price_ being
a particularly effective method (e.g., Desaulniers et al., 2006). Second,
there is also by now a substantial literature on _stochastic_ MIPs (e.g.,
Küçükyavuz and Sen, 2017). Third, considerable effort has been made to extend
the MIP algorithmic technology to cope with nonlinearities, leading to the
blossoming field of _mixed-integer nonlinear programming_ or MINLP (e.g., Lee
and Leyffer, 2012). Particularly effective algorithms and software packages
are now available for _convex_ MINLP (e.g., Kronqvist et al., 2019), and one
of its important special cases, mixed-integer second order cone programming
(e.g., Benson and Sağlam, 2013).

###  2.16 Nonlinear programming272727This subsection was written by E. Alper
Yıldırım.

Nonlinear programming is a generalisation of linear programming (§2.14), in
which the objective function or the constraints can be given by general
nonlinear functions. Mathematically, a nonlinear programming problem is
represented as

$\textrm{(P)}\quad\min\\{f(x):x\in\mathcal{S}\\}.$

Here, $\mathcal{S}=\\{x\in\mathbb{R}^{n}:g_{i}(x)\leq 0,\quad i=1,\ldots,m\\}$
denotes the feasible region, where,
$g_{i}:\mathbb{R}^{n}\to\mathbb{R},~{}i=1,\ldots,m$, and
$f:\mathbb{R}^{n}\to\mathbb{R}$ denotes the objective function.

In comparison with linear programming, nonlinear programming problems have
much more expressive power. As such, nonlinear programming problems naturally
arise in almost every setting, ranging from investment planning to machine
learning; from engineering to medicine; and from energy to sustainability
(§3.14; §3.5; §3.19; §3.9; §3.11; §3.13).

In this subsection, we will give a brief overview of theory and algorithms.
While we will not cover the modelling aspect, we will mention some classes of
optimisation problems with desirable properties, which should imply that using
an optimisation model from such classes would significantly increase the
likelihood of solving it.

The difficulty of the generic optimization problem (P) is largely determined
by the properties of the objective function $f:\mathbb{R}^{n}\to\mathbb{R}$
and of the functions $g_{i},~{}i=1,\ldots,m$ that define the feasible region
$\mathcal{S}\subseteq\mathbb{R}^{n}$. Generally speaking, increasingly more
restrictive assumptions on $f$ and on $g_{i},~{}i=1,\ldots,m$ give rise to
increasingly more structured optimisation problems with stronger and more
desirable properties. For instance, the special case in which each of $f$ and
$g_{i},~{}i=1,\ldots,m$ is a linear function, referred to as _linear
programming_ (§2.14), is arguably the most structured class of optimization
problems with very appealing theoretical properties, which lay the groundwork
for several effective solution methods such as the simplex method (see, e.g.,
Dantzig, 1990) and interior-point methods (Karmarkar, 1984; Wright, 1997; Ye,
1997). In contrast, general nonlinear programming problems usually enjoy fewer
desirable properties.

The class of _convex optimisation problems_ is comprised of optimisation
problems in which each of $f:\mathbb{R}^{n}\to\mathbb{R}$ and
$g_{i},~{}i=1,\ldots,m$ is a convex function, which implies that
$\mathcal{S}\subseteq\mathbb{R}^{n}$ is a convex set, and includes linear
programming as a special case. Any optimisation problem that does not belong
to this class is a _nonconvex optimisation problem_. On the other hand, (P) is
called an _unconstrained optimisation problem_ if
$\mathcal{S}=\mathbb{R}^{n}$, and a _constrained optimisation problem_
otherwise.

A useful notion in nonlinear programming is that of _local optimality_. A
point $\hat{x}\in\mathbb{R}^{n}$ is said to be a _local minimiser_ of (P) if
there exists an open ball $\mathcal{B}\subset\mathbb{R}^{n}$ of positive
radius centred at $\hat{x}$ such that $\hat{x}$ is a minimiser of $f$ over the
potentially smaller feasible region $\mathcal{B}\cap\mathcal{S}$. In contrast,
$\hat{x}$ is a global minimiser of (P) if $\hat{x}$ is a minimiser of $f$ over
the entire feasible region $\mathcal{S}$. Note that a global minimiser is also
a local minimiser.

We next briefly give an overview of optimality conditions for each
aforementioned class of optimisation problems. We start with unconstrained
optimisation problems in the one-dimensional setting (i.e., $n=1$). If
$\hat{x}\in\mathbb{R}$ is a local minimiser of (P), then $f$ should be neither
decreasing nor increasing at $\hat{x}$. Assuming that $f$ is a continuously
differentiable function, we therefore obtain $f^{\prime}(\hat{x})=0$. This
geometric interpretation carries over to the higher-dimensional setting (i.e.,
$n\geq 2$) by simply viewing a mutivariate function as a collection of one-
dimensional functions along _feasible directions_ at each
$\hat{x}\in\mathbb{R}^{n}$, i.e., directions along which one can move starting
from $\hat{x}\in\mathbb{R}^{n}$ and still remain in the feasible region. In
the unconstrained case, every direction $d\in\mathbb{R}^{n}$ is a feasible
direction at every $x\in\mathbb{R}^{n}$. Using the result from the one-
dimensional case, if $\hat{x}\in\mathbb{R}^{n}$ is a local minimiser of (P),
then the partial derivatives of $f$ with respect to each variable should be
zero, or equivalently, that $\nabla f(\hat{x})=0\in\mathbb{R}^{n}$, where
$\nabla f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the gradient of $f$. Such a
point is called a _stationary point_.

For the special case of convex unconstrained optimisation problems, the
convexity of the objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ implies
that the aforementioned necessary conditions are also sufficient, i.e., a
point is a local minimiser if and only if it is a stationary point.
Furthermore, for convex functions, every local minimiser is, in fact, a global
minimiser. Therefore, we obtain the equivalence between global minimisers and
stationary points. On the other hand, for a nonconvex optimisation problem,
there may be stationary points that may not correspond to a local minimiser of
(P) (e.g., if $f(x)=x^{3}$, then $\hat{x}=0$ is a stationary point but not a
local minimiser). As illustrated by this example, the complete
characterisation of global optimality does not carry over from convex
optimisation to nonconvex optimisation, even in the unconstrained setting.

For constrained optimisation problems, we first consider the convex
optimisation case. By the convexity of the feasible region
$\mathcal{S}\subseteq\mathbb{R}^{n}$, for any $\hat{x}\in\mathcal{S}$, the set
of all feasible directions is given by $\tilde{x}-\hat{x}\in\mathbb{R}^{n}$,
where $\tilde{x}\in\mathcal{S}$. Arguing similarly to the unconstrained case
and using the convexity of $f$, a point $\hat{x}\in\mathcal{S}$ is a global
minimiser of (P) if and only if $f$ does not decrease along any feasible
direction, i.e., if and only if $\nabla f(\hat{x})^{T}(\tilde{x}-\hat{x})\geq
0$ for all $\tilde{x}\in\mathcal{S}$. Therefore, as in the unconstrained case,
we once again have the equivalence between local and global minimisers.

Next, consider a nonconvex constrained optimisation problem. If the feasible
region $\mathcal{S}\subseteq\mathbb{R}^{n}$ is a convex set but $f$ is a
nonconvex function, a similar argument as in the convex case gives rise to the
following necessary condition: If $\hat{x}\in\mathcal{S}$ is a local minimiser
of (P), then $\nabla f(\hat{x})^{T}(\tilde{x}-\hat{x})\geq 0$ for all
$\tilde{x}\in\mathcal{S}$. As in the unconstrained case, simple examples show
that this condition is no longer sufficient for local optimality. If, on the
other hand, $\mathcal{S}\subseteq\mathbb{R}^{n}$ is a nonconvex set, then we
instead rely on a more general notion of _tangent directions_ to the feasible
region $\mathcal{S}$ at $\hat{x}$. Therefore, if $\hat{x}\in\mathcal{S}$ is a
local minimiser, then $\nabla f(\hat{x})^{T}d\geq 0$ for every tangent
direction $d\in\mathbb{R}^{n}$ to the feasible region $\mathcal{S}$ at
$\hat{x}$. In general, the set of such tangent directions may not be easy to
characterise. Under certain additional assumptions about the geometry of the
feasible region $\mathcal{S}\subseteq\mathbb{R}^{n}$, referred to as
_constraint qualifications_ (Bazaraa et al., 2005, Chapter 5), explicit
necessary optimality conditions can be derived.

Having reviewed optimality conditions, we finally give a brief overview of
methods for solving optimisation problems. Nonlinear optimisation algorithms
are generally _iterative_ in nature, i.e., they generate a sequence of points
$x_{k}\in\mathbb{R}^{n},~{}k=1,2,\ldots$ that satisfies certain properties.
For instance, the sequence may either converge to a local or global minimiser
of an optimisation problem, or may simply have a limit point that satisfies
the necessary conditions for local optimality. As illustrated by the
discussion on optimality conditions, one can establish considerably weaker
properties for nonconvex optimisation problems in comparison with convex
optimisation problems. In fact, most classes of nonconvex optimisation
problems are provably difficult in a formal complexity sense, even when
restricted to minimising a quadratic function over a polyhedron (Murty and
Kabadi, 1987; Pardalos and Vavasis, 1991). As such, it would not be reasonable
to expect an algorithm to solve every optimisation problem to global
optimality in a reasonable amount of time.

Therefore, different performance metrics are employed for assessing algorithms
for different classes of optimisation problems. While, for convex optimisation
problems, one usually expects a “good” algorithm to compute a global optimal
solution, an algorithm for a nonconvex optimisation problem could be deemed
“effective” if it always converges to a local (rather than a global) optimal
solution.

In the unconstrained case, given an iterate $x_{k}\in\mathbb{R}^{n}$, the main
idea is to identify a feasible direction $d\in\mathbb{R}^{n}$ along which the
objective function will decrease. Such a direction $d\in\mathbb{R}^{n}$,
called a _descent direction_ , would necessarily satisfy $\nabla
f(x_{k})^{T}d<0$. Then, a step size in this direction is determined according
to certain criteria that would guarantee a decrease in the objective function.
Therefore, this family if algorithms is referred to as _gradient descent
methods_ and includes _steepest descent_ as a special case (i.e., the case
where $d=-\nabla f(x_{k})$). Under mild assumptions, this class of algorithms
converges to a stationary point of $f$. Recall that such a point is a global
minimiser if $f$ is a convex function. Other methods in this class are Newton
methods, conjugate gradient methods, and quasi-Newton methods, each of which
generates iterates that converge to a stationary point under appropriate
assumptions.

Considering the constrained case, while general convex optimisation problems
do not retain all desirable properties of the simpler class of linear
programming problems, they still have a sufficiently rich structure that pave
the way for provably efficient solution algorithms. In fact, every convex
optimisation problem, in theory, can be solved to global optimality by the
ellipsoid method (Yudin and Nemirovskii, 1976; Shor, 1977) or by interior-
point methods in polynomial time (Nesterov and Nemirovskii, 1994).
Furthermore, a variety of highly effective commercial and non-commercial
solvers are available for solving several classes of convex optimisation
problems such as linear programming, second-order cone programming, and
semidefinite programming that frequently arises in applications (see, e.g.,
https://neos-server.org/neos/solvers/index.html).

For the nonconvex constrained case, one approach is based on approximating a
constrained optimisation problem by a sequence of unconstrained optimisation
problems by either using a penalty function, based on penalising violation of
constraints (_penalty methods_), or using a barrier function, based on
preventing the violation of constraints by keeping the iterates strictly in
the relative interior of the feasible region
$\mathcal{S}\subseteq\mathbb{R}^{n}$ (_barrier methods_). Other methods
include Augmented Lagrangian methods, based on combining Lagrangian relaxation
with penalty methods, and Sequential Quadratic Programming methods, based on
approximating the optimisation problem by a quadratic programming problem.

Finally, various real-life applications in machine learning and data science
give rise to very large-scale problems that are beyond the capability of
current solvers and computing platforms. For such problems, there exist a
variety of heuristic optimisation methods that can be employed to find a good
solution in a reasonable amount of time (§2.13). However, in contrast with
exact methods, such methods usually do not provide any guarantees on the
quality of the solution.

Nonlinear optimisation is a very active area of research. The reader is
referred to excellent textbooks for further information (e.g., Fiacco and
McCormick, 1968; Mangasarian, 1994; Bazaraa et al., 2005; Nocedal and Wright,
2006; Bertsekas, 2016; Luenberger and Ye, 2016).

###  2.17 Queueing282828This subsection was written by Hayriye Ayhan and Tuğçe
Işık.

Queueing systems arise in many real life applications including production,
service systems, finance, logistics and transportation. As mentioned in
Stidham (2002), many queueing models have been studied even before the
introduction of Operational Research in the 1950s. We provide a brief overview
of methodologies used in queueing systems analysis. We start with exact
methods and then continue with approximations and asymptotic analysis.

The classical analysis of queueing systems involved modelling the single stage
Markovian queues as birth-death processes and computing their steady state
performances using Markov chain theory. These earlier theoretical
contributions were initially summarised in Feller’s two volume books Feller
(1957a, b), then in classical textbooks such as Cooper (1972), Gross and
Harris (1974) and Kleinrock’s two volumes Kleinrock (1975a, b), and more
recently in Gautam (2012), Harchol-Balter (2013), and many other books. Takács
(1962) focused on using transforms and generating functions for steady state
and transient behaviour of queueing systems. In the early days, transforms and
generating functions were considered to be exact expressions but one had to
invert these generating functions in order to obtain the actual performance
measures which is in general difficult. Marcel Neuts was the first one who
approached this inversion problem algorithmically. In his 1981 book, Neuts
(1981) focused on queues that generalise the G/M/1 structure, whereas in his
second book, Neuts (1989) generalised the structure of the M/G/1 queue. The
main idea in these books is to approximate the non-exponential distributions
with a phase type distribution (convolution and mixture of exponentials) which
yields a continuous time Markov chain model for the original system that could
be analysed, at least numerically. This line of research resulted in many
contributions on the so-called matrix-geometric methods (see also Latouche and
Ramaswami, 1999). Arguably the most well known result in queueing theory is
Little’s law ($L=\lambda W$ or its generalisation $H=\lambda G$) which
provides a relationship between the mean steady state number of customers and
the mean sojourn time in a system. For a thorough survey of the Little’s
result and its extensions, the reader is referred to Whitt (1991). There are
numerous proofs of Little’s law but El-Taha and Stidham (1999) provide an
elegant sample path proof. On the other hand, Bertsimas and Nakazato (1995)
relate the steady-state distribution of the number in the system (or in the
queue) to the steady state distribution of the time spent in the system (or in
the queue) in a queueing system under FIFO (First In First Out).

While there has been a lot of interest in stationary queues, Massey’s 1981
dissertation Massey (1981) drew attention to the analysis of non-stationary
queues (i.e., queues with time dependent arrival and service processes).
Massey’s dissertation started with the analysis of M(t)/M(t)/1 queue and then
extended to other non-stationary Markovian systems. Many subsequent papers,
such as Massey and Whitt (1998), focused on queueing models with time-
dependent arrival rates, especially infinite-server “offered-load” models
which describe the load that would be on the system if there were no limit to
the available resources. The main idea of these papers is to provide
algorithms (approximations) to solve the Poisson equation. On the other hand,
Bertsimas and Mourtzinou (1997) derived a set of transient distributional laws
that relate the number of customers in the system (queue) at time t to the
system (waiting) time of a customer that arrived to the system (queue) at time
$t$.

Networks of queues have been of interest to researchers since 1950s. Jackson
(1957) was the first one to observe that joint steady state distribution of
the number of customers at the nodes of a network of Markovian queues with
single server (at each node) is the product of individual distribution of
M/M/1 queues. Jackson (1963) generalised this result to networks of queues
with multiple servers at the nodes. Gordon and Newell (1967) discovered that
the stationary distribution again has a product form in closed Markovian
networks but in this case a normalisation constant is required. Baskett et al.
(1975) proved that the product form is insensitive to the service time
distribution if the service discipline satisfies certain assumptions. This and
other insensitivity results in networks were also considered by Kelly (1979)
and Serfozo (1999) which also has results on other networks such as those with
blocking and rerouting. Daduna (2001) focused on obtaining explicit
expressions for the steady behaviour of discrete time queueing networks and
gave a moderately positive answer to the question of whether there can be a
product form calculus in discrete time. In recent years, a number of models
involving different compatibilities between jobs and servers in queueing
systems, or between agents and resources in matching systems, have been
studied, and, under Markovian assumptions and appropriate stability
conditions, the stationary distributions were again shown to have product
forms (see Gardner and Righter, 2020, and the references therein). Baccelli et
al. (1992) modelled a class of networks using the so-called (max,+) linear
systems. In their pioneering work, using (max,+) algebra techniques, Baccelli
and Schmidt (1996) derived Taylor series expansions for the mean waiting times
in Poisson driven queueing networks that belong to the class of (max,+) linear
systems. Even though these expansions are sometimes referred to as light
traffic approximations, in some cases all coefficients of the series expansion
can be computed yielding an exact expression. These results were generalised
to transient performance measures by Baccelli et al. (1997) and joint
characteristics by Ayhan and Baccelli (2001).

Exact analysis of general queueing systems is often challenging, making the
characterisation of performance measures difficult. Thus, asymptotic analyses
are commonly carried out via various approximation methods. We next provide an
overview of such methods.

Many of the earlier works on the asymptotic analysis of queueing systems
focused on heavy traffic and many server approximations for single stage
queues. In his pioneering work, Kingman (1961, 1962, 1965) have asymptotically
characterised the waiting time distribution for the single server queue with
general interarrival and service time distributions under heavy traffic
conditions (i.e., when the traffic intensity $\rho\to 1$). Several others have
developed heavy traffic approximations for the G/G/s queue, where a sequence
of systems with fixed number of servers and traffic intensities
$\\{\rho_{n}\\}$ approaching one are considered (see, for example,
Köllerström, 1974). In these approximations, the sequence of normalised (i.e.,
scaled) queue length processes converge to a reflected Brownian motion with
negative drift (see Whitt, 2002), and the associated sequence of scaled
stationary queue-length distributions (i.e., the stationary distribution of
the limiting diffusion process) converges to an exponential distribution. We
refer the reader to Harrison (1985) for a detailed technical treatment of
heavy traffic limits and diffusion approximations. Asymptotic analysis was
also considered for multi-class and multi-stage queueing networks. Defining
the stability region of these networks using fluid limit analysis was
considered in Chen (1995); Dai (1995); Dai and Meyn (1995). Many of the works
that considered heavy traffic analysis of multi-class queueing networks focus
on achieving the so-called state space collapse. Bramson (1998) demonstrated
the state space collapse for first-in first-out queueing networks of Kelly
type and head-of-the-line proportional processor sharing queueing networks.
His framework has been used to prove state space collapse results in several
other works including Stolyar (2004) and Mandelbaum and Stolyar (2004). For a
more comprehensive review of heavy traffic analysis of multi-class queueing
networks, we refer the reader to Chen and Yao (2001).

Many-server approximations were also considered for asymptotic analyses of
queueing systems. In these approximations, the traffic intensity can be kept
constant while letting the arrival rate and the number of servers go to
infinity. Iglehart (1965) showed that the resulting sequence of normalised
queue length processes converges to an Ornstein-Uhlenbeck process in the many
server setting when the service time distributions are exponential. Later on,
Whitt (1982) generalised this result for systems with non-exponential service
times. For a more comprehensive overview of results in this area, see Whitt
(2002). In their seminal work, Halfin and Whitt (1981) defined the so called
Halfin-Whitt regime for the GI/M/s queue where the traffic intensities
converge to one from below, the number of servers and arrival rates tend to
infinity, but steady-state probability that all servers are busy remains
fixed. They showed that under the appropriate scaling, the queue length
processes converge to a diffusion process. In the past decades, many other
asymptotic results have been obtained for many server queues in the Halfin-
Whitt regime. Reed (2009) studied the G/GI/s queue and obtained fluid and
diffusion limit results for the queue length process. We refer the reader to
van Leeuwaarden et al. (2019) for a further review of the various asymptotic
results obtained in the Halfin-Whitt regime.

Although heavy traffic approximations for queues have been popular in recent
decades, light traffic (as the traffic intensity $\rho\rightarrow 0$) and
interpolation approximations have also been developed. Bloomfield and Cox
(1972) developed light traffic approximations for a single server queue.
Burman and Smith (1983) developed approximations for the expected delay in
M/G/s queue both for heavy and light traffic, and showed that as traffic
intensity goes to zero, probability of delay depends only on mean service time
distributions. Daley and Rolski (1992) used light traffic approximations to
study the limiting properties of the waiting time in many-server queues. Light
traffic approximations have also been used to study the limiting processes in
queueing networks (see, for example, Simon, 1992).

As mentioned earlier, approximation methods were commonly used in the
asymptotic analysis of time-varying (i.e., non-stationary) queues. Mandelbaum
et al. (1999) developed a fluid approximation for the queue length process in
time-varying multiserver queue with abandonments and retrials. Pang and Whitt
(2010) have developed heavy traffic approximations for infinite server queues
with time-varying arrivals. The reader is referred to Whitt (2018) for a
recent review of the literature on non-stationary queues.

Due to the interest in communication/telecommunication systems, in late 1990s
and early 2000s, there was a lot of research on queues with heavy tailed
interarrival and/or service times. Intuitively, heavy tailed distributions
decay slower than an exponential distribution (see Resnick, 2007, for a
thorough discussion). Boxma and Cohen (2000) provided an overview of results
for single service queues with heavy tailed interarrival and/or service time
distributions. Baccelli et al. (1999) and Ayhan et al. (2004) showed that the
asymptotics of response time was dominated by the station with the heavy
tailed service time in a class of open and closed networks, respectively. Foss
and Korshunov (2012) developed upper and lower bounds on the tail distribution
of the stationary waiting time in the GI/GI/s queue with heavy tailed service
times.

###  2.18 Risk analysis292929This subsection was written by Louis Anthony Cox,
Jr.

Risk analysis is a discipline that seeks to inform people about what might
happen and how to reduce the probability and severity of undesired outcomes.
It draws on decision analysis, game theory, and other areas of Operational
Research but is distinguished from them by the questions it asks, the
frameworks it provides for answering them, and the uses to which its answers
are put (Aven, 2020; Greenberg et al., 2020). Where decision analysis focuses
on principles for identifying logically coherent choices that make preferred
outcomes more likely based on a decision maker’s beliefs and value trade-offs,
risk analysis seeks to inform analytic-deliberative decision-making by
multiple stakeholders –- possibly with conflicting worldviews, values, and
beliefs –- for managing critically important matters ranging from the safe
operation of nuclear power plants to priority-setting for public and
occupational health and safety measures.

Risk analysis is often subdivided into risk perception, risk assessment, risk
communication, risk management, and risk governance and policy-making
(Greenberg and Cox, 2021). The following sections describe these components.

#### Risk Perception

Public concerns and political appetite to address them are shaped by perceived
risks, whether or not they are accurate. Several frameworks have been
developed to help understand the technical, psychological, and social drivers
of risk perceptions (Siegrist and Árvai, 2020). The psychometric paradigm
(Slovic, 2000) explains many aspects of risk perceptions in terms of a few
underlying factors such as as dread risk (associated with a lack of control,
dreaded consequences, catastrophic potential, inequity in the distribution of
risks, risks increasing over time, and fatal consequences) and unknown risk
(associated with unobservability, novelty, unknown exposure, being unknown to
science, and delayed consequences). The cognitive heuristics and biases
literature positions risk perceptions within a “dual process” framework in
which rapid emotional evaluations (“System 1”) can be modified by slower, more
effortful cognition (“System 2”) (Kahneman, 2011; Skagerlund et al., 2020).
The cultural theory of risk (Douglas and Wildavsky, 1983; McEvoy et al., 2017;
Bi et al., 2021) posits that individual perceptions of risk are shaped by
social and ideological processes that emphasise or suppress perceptions of
risks depending on the respondent’s values and preferred form of social order.
The social amplification of risk framework (SARF) (Kasperson et al., 2022)
describes the social amplification or attenuation of perceived risks as risk
information is communicated among people with different worldviews.

Major lessons from the study of risk perception are that experts and members
of the public often view risks quite differently; that experts often focus on
the probability and consequence severity dimensions of risk while members of
the public consider many other aspects; that most people tend to overestimate
the frequencies of rare but vivid events (e.g., terrorist attacks, murders)
and underestimate the frequencies of common but familiar ones (e.g., car
accidents, heart attack fatalities); and that risk perceptions of both experts
and lay people are predictably shaped and distorted by cognitive heuristics
and biases and are amplified or attenuated by media reports and other
communications in ways that reflect the recipients’ worldviews. System 1 tends
to be innumerate, responding emotionally to possibilities and categories of
harm while underweighting or ignoring relevant frequencies and magnitudes.
System 2 often fails to sufficiently adjust or correct the promptings of
System 1 leading to decisions with predictable regrets. These findings help to
explain why expert and actuarial assessments of risk often differ from lay
perceptions of risk. In a democratic society, perceptions affect decisions. A
major challenge for risk analysis is to assess and communicate risks to help
inform and improve collective decisions in ways that understand and respect
the realities of risk perception.

#### Risk Assessment

Risk assessment addresses how large and uncertain risks are. It begins with
qualitative questions about what might go wrong and proceeds to quantitative
assessments of how likely adverse events are to occur and what their possible
consequences and their probabilities would be (Kaplan and Garrick, 1981).
Probabilistic risk assessment (PRA) and quantitative risk assessment (QRA)
methods apply probability models and statistical methods to data and modelling
assumptions to quantify or bound the predicted frequencies and severities of
losses and to estimate how their joint probability distribution would be
changed by different risk management policies or interventions. Quantitative
measures of risk can be derived from the full probability distribution or
stochastic process descriptions of uncertain outcomes (Smidts, 1997),
including dynamic coherent risk measures used in financial risk analysis
(Bielecki et al., 2017). Stochastic process models of the occurrence
frequencies of adverse events (such as accidents at power plants or tornadoes
in cities) and the probability distribution of losses for each event can also
be used to estimate entire cumulative probability distributions for losses
over a stated time interval for different scenarios or sets of assumptions
(Kaplan and Garrick, 1981, see also §2.10 and §3.19).

In the past decade, PRA techniques such as causal Bayesian networks (BNs),
dynamic Bayesian networks (DBNs), and related probabilistic graphical models
have increasingly been used to predict the probabilistic effects caused by
interventions in engineering systems (Ruiz-Tagle et al., 2022) and public
health applications (Butcher et al., 2021). They have largely supplanted older
and less general PRA techniques such as fault tree analysis and Markov
decision processes (Hanea et al., 2022; Cox et al., 2018). Together with
discrete-event stochastic simulation models and continuous (systems dynamics)
simulation, they provide constructive methods to predict how risk management
interventions would change the probabilities of outcomes over time. This
information can be used for simulation-optimisation of risk management
decisions (Better and Glover, 2011).

PRA techniques have been extended to address adversarial risks in which
intelligent adaptive adversaries rather than chance events threaten the safety
and values that a risk manager seeks to protect (Banks et al., 2022); and
unknown risks (or risks under radical uncertainty, sometimes called Knightian
uncertainty by economists), in which relevant probabilities are unknown, e.g.,
by using uncertainty sets that replace precise probability distributions by
(usually convex) sets of possible probability distributions or by scenarios of
possibilities that are not necessarily exhaustive (Gilboa et al., 2017).
Recent artificial intelligence and machine learning (AI/ML) methods are now
being applied to natural hazards and disasters (Guikema, 2020), cybersecurity
(Nifakos et al., 2021), power markets (Marcjasz et al., 2022), and financial
portfolio risk management problems where new, changing, and unknown conditions
make it necessary to learn effective risk prediction and management decision
rules from data and experience without the guidance of well-validated PRA
models (Cox, 2020).

#### Risk Management, Governance, Communication, and Risk-Cost-Benefit
Analysis

Given public perceptions and technical estimates of risks, what should be done
about them? Who should decide, and how? Managing risks to human health,
safety, or the environment often involves “wicked” decision problems and
“deep” uncertainties, meaning that there are no clear, widely agreed-to
definitions of the decision problem and solutions to it (Lempert and Turner,
2021). Although multiobjective and risk-sensitive or risk-constrained
optimisation problems can be formulated for some risk management problems,
such as routing hazardous cargo, for many wicked risk management problems,
relevant decision variables, constraints, possible outcomes, and objective
functions may be unknown or not widely agreed to. Risk management in such
challenging cases usually involves issues of causation (what can be done and
how much difference in outcome probabilities would different feasible choices
make?), collective choice (how should the disparate perceptions and
preferences of individuals be resolved or aggregated for purposes of
collective decision-making?) and risk governance (who should be responsible
for making, implementing, obeying, enforcing, and revising risk management
decisions; how should stakeholders participate in risk management decisions;
what institutions and processes should guide, restrain, and integrate
collective risk management; and how should conflicts be resolved and
collective decisions be made at individual, organisational, community, local,
national, and international levels?) (Klinke and Renn, 2021).

The most immediate decisions for risk managers responding to potential or
actual crises are often about risk communication. For example, if a pandemic
or natural disaster such as a tsunami, volcanic eruption or hurricane, seems
possible but not necessarily imminent or certain, then what should scientists
and government officials tell policy-makers and the public about the uncertain
risks? Who should say what to whom, and how soon? Different risk communication
goals such as informing and empowering individual and community decisions,
persuading individuals to change their behaviours, instructing citizens what
to do, and informing or shaping policy deliberations and decisions, require
different communication approaches. Risk communication frameworks for
addressing these challenges overlap with risk perception frameworks but also
emphasise the roles of trust in information sources and of outrage in
mobilising public engagement and changing behaviours (Malecki et al., 2020).

Risk-cost-benefit analysis provides a simple-sounding approach to collective
risk management: take risk management actions to maximise expected social
utility or net societal benefits (expressed as expected net present value,
possibly with a risk-adjusted discount rate; Eliashberg and Winkler, 1981;
Hammond, 1992). However, mathematical impossibility theorems have shown that
when different people have sufficiently different beliefs and preferences,
there may be no coherent way to aggregate them to make collective decisions
that respect normative principles such as Pareto efficiency (Nehring, 2007).
Trade-offs between measures of accuracy and fairness have recently been
identified for machine-learning algorithms used in risk assessment in areas
such as mortgage lending and criminal justice (Corbett-Davies et al., 2017).
Societal risk management is now often viewed less as a top-down or centralised
decision and control process in which experts provide estimated probabilities
and social utilities or net benefits to use in risk-cost-benefit or social
utility maximisation calculations than as a participatory democratic
deliberative process (Rad and Roy, 2021). Experts in risk analysis can provide
useful technical information in this process but should not dominate it
(Pellizzoni and Ungaro, 2000; Greenberg and Cox, 2021; Klinke and Renn, 2021).

Risk analysis poses intellectual, technical, and practical implementation
challenges that are likely to engage and challenge Operational Research and
risk analysis professionals for the foreseeable future. A more detailed review
of the accomplishments, current state, and remaining challenges for much of
the field of risk analysis can be found in the 40th Anniversary Special Issue
of Risk Analysis (Greenberg et al., 2020), in specialised books such as Aven
(2015), and in online resources303030For example, https://www.sra.org/risk-
analysis-introduction/.

###  2.19 Simulation313131This subsection was written by Christine S.M.
Currie.

A simulation aims to reproduce the important behaviour of a real system. Our
focus here is on the use of computer simulation models within operational
research (OR), whilst acknowledging that the field is much wider and ranges
from computer models of sub-atomic particles to simulations involving real
human actors, particularly prevalent in medicine and health sciences. Three
main flavours of simulation are used in this context: discrete event
simulation, agent-based modelling and system dynamics. After discussing the
uses of simulation, we continue this subsection by introducing these three
main flavours before going on to discuss four important areas in simulation
research: conceptual modelling, input modelling and parameterisation,
simulation optimisation, and finally the newer area of data-driven simulation,
linking to Industry 4.0 and digital twins. A selection of open source tools
for simulation are given in §3.18.

A simulation model, built on a computer, has a number of potential functions.
Principally it is used for experimentation because testing out new settings or
ways of working on a simulation results in fewer negative implications than
experimenting with the real system. This can allow simulation to be used for
optimisation of complex stochastic systems and there has been considerable
research in this area in recent years, as we discuss below. Simulation can
also be used for predicting future behaviour, and the COVID-19 pandemic
showcased the predictive power of simulation modelling in a very high-profile
situation (e.g., the agent-based model used to advise the UK government and
described in Ferguson et al., 2020). The process of building a simulation
model results in a better understanding of the real system because of the need
to identify and model the important relationships between different entities.
Running the model can also help with estimating the sensitivity of model
outputs to system parameters. Another use of simulation models is for
training. Within the OR context, this most often takes the form of strategic
game-playing (e.g., the beer game developed by MIT and described in Sterman,
1989) to practice decision-making under different scenarios in a safe
environment.

Discrete event simulation (DES) is typically used to model systems in which
entities move through a set of activities. Where these activities require
resources, entities will queue until the resource becomes available. Such
simulations are described as discrete event because the system state only
varies at discrete time points, known as events. For example, a DES model
might be used to describe a manufacturing line and in this case the events
could include an item starting or finishing processing by a machine on the
line. Usually in DES, the simulation clock will jump from one event to the
next rather than moving in equal time steps.

System dynamics (SD; §2.22) was first developed in the 1950s by Jay Forrester
to help with the understanding of industrial problems. SD models deal with
stocks and flows, where the dynamics are dictated by a set of connected
differential equations. Describing a system using an SD model can help with
detecting feedback loops and delay effects and SD modelling is useful for
strategic decision-making.

Agent based modelling (ABM) describes the behaviour of individual entities or
agents within a population. As Macal (2016) states, one of the key differences
between ABM and both DES and SD is that it takes an agent perspective of a
system. Each agent in the simulation will follow its own set of rules
dictating its behaviour and how it interacts with the environment and other
agents. Agent behaviour is typically stochastic, allowing natural variability
to be included in the model. An ABM can be used as a bottom up approach to
determine emergent behaviour where individual actions lead to a system level
response.

Regardless of the simulation modelling technique, the first part of any
simulation project is to gain an understanding of the system being modelled,
the objectives of the work, and the key components that should be included,
referred to as conceptual modelling. There is some discussion of the exact
definition of conceptual modelling in Robinson (2008) but some key points are
made to support the process. First, the conceptual model can be thought of as
separate from the final computer model that is built and serves as an
abstraction of the real system that describes what is going to be modelled and
gives an indication of how that might be done. Second, the development of the
conceptual model requires input from the modeller and the system owners.
Third, a conceptual model does not remain constant through a simulation
project but is revisited and adapted as the project continues. Recent research
in conceptual modelling is reviewed by Robinson (2020) and has focused on
designing modelling frameworks. The work described tends to be related to DES
models but the core principles can also apply to the building of both ABS and
SD models.

Like any other model, the utility of a simulation is very much dependent on
its inputs: the garbage-in-garbage-out principle holds true here. Setting up
the probability distributions used for inputs of a stochastic simulation model
or parameterisation of a deterministic simulation model is referred to as
input modelling and is typically achieved through fitting statistical models
to available data and eliciting expert opinion. When estimating the inputs for
a simulation model from data there is some uncertainty in their true values.
With a different set of data, the estimates of the inputs would likely be
different. Any uncertainties setting the model inputs will propagate through
to the model outputs, resulting in input uncertainty. This is influenced both
by the accuracy of the estimates of the inputs and the sensitivity of the
model output to that particular distribution or parameter. Corlu et al. (2020)
provide a review of the current state of the art in input uncertainty research
for simulation, while Song et al. (2014) provide practical suggestions on how
to estimate the impact of unput uncertainty on the output results.

Often simulation models are used to experiment with different system set-ups.
Simulation optimisation, sometimes referred to as optimisation via simulation
describes the use of a simulation model to find the optimal value for one or
more decision variables. Typically it is used in the design of stochastic
systems that are too complex to be effectively described by an analytical
model. Practical examples of problems that can be solved using simulation
optimisation include finding the optimal number and configuration of beds in a
hospital ward; determining the appropriate number of repair staff on a
production line; choosing between a selection of different configurations for
a system.

The problem can be represented mathematically as

$\min g(\mathbf{x}),\mbox{ }\mathbf{x}\in\mathbf{\Theta},$

where the function we are optimising $g(\mathbf{x})$ is generally the expected
value of the output of a stochastic simulation model,
$g(\mathbf{x})=E[Y(\mathbf{x},\xi)]$; $\mathbf{x}$ is a vector of decision
variables; $\mathbf{\Theta}$ denotes the feasible region for $\mathbf{x}$; and
$\xi$ indicates the randomness inherent in the model. The majority of research
in simulation optimisation aims to improve the efficiency of the optimisation
algorithms by reducing the number of simulation replications needed to
estimate the optimal values of $\mathbf{x}$. Where a complex and slow-running
simulation model is used to generate the $Y(\mathbf{x},\xi)$ this efficiency
is particularly important. Hong and Nelson (2009) classify simulation
optimisation problems into three categories:

  1. 1.

Small number of solutions: $\mathbf{\Theta}$ contains only a small number of
solutions and the decision variable $\mathbf{x}$ might define a particular
system configuration. In this case the problem is one of ranking and
selection.

  2. 2.

Decision variables are continuous: $\mathbf{\Theta}$ is a convex subset of the
set of $d$-dimensional real numbers and the problem is continuous optimisation
via simulation.

  3. 3.

Decision variables $\mathbf{x}$ are discrete and ordered: $\mathbf{\Theta}$ is
a subset of the set of $d$-dimensional integers and the problem is discrete
optimisation via simulation.

A set of algorithms exists for solving each category of problem. There has
also been significant work on multi-objective optimisation via simulation; for
example, see Hunter et al. (2019) for a detailed description of the problem
and different solution approaches.

In recent years, sensing has become more widespread and the transfer of data
from physical systems to control systems is now happening in close to real
time. This has allowed modellers to design simulation models that are
automatically fed data from the real system allowing them to either predict
the future (e.g., prediction of emergency departments crowding Hoot et al.,
2008) or to use the simulation models to optimise system parameters as part of
a dynamic control process. Such models are sometimes referred to as a digital
twin or symbiotic simulation (see Onggo, 2019, for a definition). Xu et al.
(2016) describe how simulation can be incorporated into the Industry 4.0
framework using the example of a semiconductor fab operation. The use of
simulation in dynamic control is still in its infancy and requires fast and
reliable simulation optimisation algorithms as well as mechanisms for enabling
the simulation model to evolve based on new input data. Industry 5.0 is
intended to complement Industry 4.0 by putting societal goals at the heart of
industrial decision-
making323232https://msu.euramet.org/current_calls/documents/EC_Industry5.0.pdf.
This is a potential growth area for simulation optimisation, particularly
multi-objective simulation optimisation (e.g., see Hunter et al., 2019) which
enables solutions to be found that support several competing objectives.

Being a huge topic with many different facets, there is no single article that
provides an overview of simulation but there are several excellent textbooks
covering simulation techniques including Law (2015) and Banks et al. (2004).
The archive of the Winter Simulation Conference333333www.wintersim.org is also
an extensive resource in understanding the state-of-the-art in the field and
its tutorial papers provide more basic tuition in the effective implementation
of simulation. Recently, the history track at the conference has also provided
an overview of the evolution of the simulation field.

###  2.20 Soft OR343434Henceforth in this subsection we just refer to PSMs to
avoid the obvious dual of ‘Hard OR’ and therefore to manufacture, or at least
to continue to propagate, an unhelpful distinction. We retain its use here for
signposting. and problem structuring methods353535This subsection was written
by Mike Yearworth and Leroy White.

Problem Structuring Methods (PSMs) are concerned with addressing problem
formulation in OR. Following definitions of Mingers and Rosenhead (2004) and
Rosenhead (1996) they consist of a set of rigorous but not mathematical
methods based on qualitative, diagrammatic modelling. They allow for a range
of distinctive stakeholder views of a problem to be expressed, explored and
accommodated. They encourage active participation of stakeholders in the
modelling process, through facilitated workshops and the cognitive
accessibility of the modelling approach. PSMs afford negotiating a joint
agenda and ownership of actions. The aim is for exploration, learning and
commitment from stakeholders, rather than optimisation or prediction. PSMs
thus are vital and constitute a significant developmental direction for OR.
See Smith and Shaw (2019) and Franco and Rouwette (2022) for recent reviews.

Understanding the contribution of PSMs to OR requires some knowledge of their
evolution. We characterise the development of the field into three phases: (i)
origins, (ii) growth (only noted here through the increased publication rate
of PSM related articles), and (iii) maturity, covering the diffusion of PSMs
to fields outside of OR, and re-integration of problem structuring into
mainstream OR. Looking at the last first, we see PSMs at an important turning
point, as recent work by Dyson et al. (2021) specifically identify the
centrality of problem structuring in the origins of OR and lead us towards the
important question of why PSMs are not seen as an essential element of every
OR engagement.

The origins of PSMs as a set of formal methods in OR arose as a consequence of
the broad critique of the process of OR in the 1970-80s; the label itself was
pioneered by Woolley and Pidd (1981). Ackoff became a trenchant critic of the
sole pursuit of objectivity and optimisation in OR describing it as an “opt-
out” (Ackoff, 1977) and set out an agenda for reconceptualising OR practice
(Ackoff, 1979a). Dando and Bennett (1981) described the situation as a
“Kuhnian crisis in management science”. In Rosenhead’s “Rational Analysis for
a Problematic World: Problem Structuring Methods for Complexity, Uncertainty
and Conflict” their prescription in OR engagements was associated with dealing
with problem contexts identified variously as wicked, messy, or swampy
(Rosenhead, 1989, pp. 3-11). These can be summarised as problem situations
that are not well-defined, involving many interested parties with different
perspectives (worldviews), where there is difficulty agreeing objectives and
the meaning of success, and that require creating agreement amongst the
parties involved for actions to be taken. The implication of the dichotomous
framing of problem contexts – i.e., wicked/tame, swamp/high ground, hard/soft,
tactical/strategic – was to set out a clear critique for the whole field of OR
and to suggest that to retain its relevance in dealing with the messiness of
real-world problems PSMs were required to bring some rigour – and indeed a
reminder of the importance – to the process of problem formulation.
Importantly, the pioneers of PSMs were concerned that traditional (i.e., ‘Hard
OR’) processes for problem formulation were practitioner-free (Checkland,
1983; Rosenhead, 1986).

The main PSMs set out by Rosenhead (1989) were Strategic Options Development
and Analysis (SODA; Eden, 1989), arising from cognitive mapping; Soft Systems
Methodology (SSM; Checkland, 1989), emerging from the failure of Hard Systems
Thinking approaches (e.g., Systems Engineering, RAND-style Systems Analysis)
when applied to messy problems; and the Strategic Choice Approach (SCA;
Friend, 1989), arising from planning. In addition, Robustness Analysis,
Metagame Analysis, and Hypergame Analysis were also included. However, setting
the boundary of PSMs has always been an open question (Mingers, 2011b). In the
main, the core methods (SODA, SSM, SCA) are seen as exemplary and provide
sufficient coherence for OR scholars and practitioners to be provided with a
clear view of a common theme.

The methodology of PSMs has long been associated with contextual matters
through Systems Thinking (§2.23 Checkland, 1983; Mingers and White, 2010),
Community OR (Johnson et al., 2018; Jones and Eden, 1981; Parry and Mingers,
1991), and large group processes (Shaw et al., 2004; White, 2002).
Methodological individualism has been addressed through Behavioural OR (§2.2;
Franco and Hämäläinen, 2016; White, 2016). There has also been long-standing
relevance to Multi-Criteria Decision Analysis (MCDA; Marttunen et al., 2017),
value-focused thinking (Keeney, 1996b), policy analysis (Eden and Ackermann,
2004), and strategy making (Ackermann and Eden, 2011; Dyson, 2000). Bridging
between PSMs and other techniques in OR has been developed as multimethodology
(Mingers and Brocklesby, 1997); for example, integration with Simulation
(§2.19; Kotiadis and Mingers, 2006; Tako and Kotiadis, 2015). Some approaches
to using the Viable Systems Model (VSM), System Dynamics (§2.22), and Decision
Analysis (§2.8) would also be considered as PSMs too (Rosenhead and Mingers,
2001, pp. 267-288), e.g., VSM (Lowe et al., 2020) and System Dynamics (Lane
and Oliva, 1998). We also see developments in Group Model Building (GMB) from
the System Dynamics community making a significant contribution to PSMs
(Andersen et al., 2007). In their growth and mature phase, applications of
SSM, SCA, and SODA have extended the reach of PSMs into, e.g., project
management (Franco et al., 2004) and environment, sustainability, and energy
policy, e.g., SCA (Fregonara et al., 2013), SODA (Hjortsø, 2004), SSM (Pahl-
Wostl, 2007), and the Drivers, Pressures, State, Impact and Response framework
(DPSIR Bell, 2012).

The state of the art and research agenda for PSMs has been the subject of
periodic reflection e.g., reviews by Rosenhead (1996) and Mingers and
Rosenhead (2004). A Special Issue of JORS in 2006 questioned where PSMs were
heading (Rosenhead, 2006) – variously argued as a “grassroots revolution”
(Westcombe et al., 2006), an appeal to common principles (Eden and Ackermann,
2006), and observations that “form and content have evolved through
interaction between the ideas and their practical use” (Checkland, 2006). A
more recent viewpoint debate “whither PSMs” again questioned their direction
of travel Harwood (2019); Lowe and Yearworth (2019).

The qualitative nature of PSM methods raises questions about evaluating both
effectiveness and value. Midgley et al. (2013), White (2006), and Franco and
Rouwette (2022) have addressed the question of effectiveness, and whilst White
goes some way towards defining the value of PSMs it is important to note the
reservations expressed by Checkland and Scholes (1990, p. 299) – that
measuring value in a unique problem context, the ‘messy’ realm of PSMs, is
unlikely to be meaningful. Tully et al. (2019) examine this conundrum in depth
from the perspective of a consulting business and make some practical
suggestions for its resolution.

Theory provides an important basis for PSM development. The range of PSM
practice reported has been explained by the constitutive rules that underpin
specific methods. First articulated by Checkland (1981, pp. 252-254),
constitutive rules are generative of method rather than prescriptive and
account for the range of practices that emerge, even when adopting a specific
methodology such as SSM i.e., adaptation is always necessary to address the
specifics of the application context. The idea was developed further by
Jackson (2003, pp. 307-311) and then by Yearworth and White (2014) into a
generic constitutive definition for PSMs. Another significant development has
been a focus on PSMs as practice and drawing on practice theories. These
theories provide a means of understanding OR practices by “zooming-in” to the
detailed, fine-grained, scale and by “zooming-out”, looking at how specific
practices affect the broader context (Ormerod et al., 2023b). Together these
theoretical strands provide sufficient basis on the one hand, to liberate PSMs
from the pigeon-hole of the dichotomous framing of their origins, and on the
other, to address OR practice as a whole and to see problem structuring as a
normal, indeed necessary, part of every OR intervention. For instance, Actor
Network Theory (ANT) provides a lens to look at the processes of
problematisation (i.e., problem formulation) in OR practice (White, 2009).
Callon (1981) draws specific attention to the “abundance of problematisations”
facing expert practitioners – that there is no single specific way of
problematising. Strands of ANT focus on the performative idiom; (Ormerod,
2014a) draws attention to the “mangle of practice” and the need for more
informative case studies of OR practice. Other theoretical underpinnings are
relevant to PSM developments e.g., PSMs as technology (Keys, 1998), Critical
Realism (Mingers, 2000), Activity Theory (White et al., 2016), and the
specific role of models as boundary objects (Franco, 2013) in facilitated
workshops (Franco and Montibeller, 2010).

From a practitioner point of view, the recent report “Reinvigorating Soft OR
for Practitioners” by Ranyard et al. (2021) to the Heads of OR and Analytics
Forum (HORAF) and the inclusion of the knowledge requirement “How to select
and apply, a range of problem structuring methods to understand complex
problems” in the Operational Research Specialist Degree Apprenticeship
specification by the Institute for Apprenticeships & Technical Education
(2021) are a welcome development.

In conclusion, for PSMs we see a return to the roots of OR as a discipline –
encompassing both practice and academic scholarship – through the centrality
of problem formulation to the process of OR (Churchman et al., 1957, p. 13)
and a reminder that the seeds of problem structuring can be seen in the work
of the ‘founders’ of OR as uncovered by Dyson et al. (2021). We have
identified a number of research gaps that indicate future research directions
for the development of PSMs. In the area of the impact of new digital
technologies, Yearworth and White (2019) propose greater use of online “same
time/different places” problem structuring workshops in order to meet the
requirements for fast meeting setup times, reducing carbon emissions, enabling
the scale-up to large group participation, and supporting new post-pandemic
working patterns. The need to address complex policy issues in the context of
wicked problems is highlighted by Howick et al. (2017) and Ferretti et al.
(2019), who argue for a re-invigorated engagement for PSM practice in policy
analysis. Finally, Ormerod (2014b), Ranyard et al. (2015) and Ormerod et al.
(2023a) remind us that we need to see a renewed practitioner-led orientation
for OR scholarship that grounds future development in solid empirical work.

###  2.21 Stochastic models363636This subsection was written by Haitao Li.

Many decision problems involve uncertainty, e.g., network design with
disruption risk, portfolio selection with uncertain return, resource planning
with unknown resource availability, crop planting with uncertain yield,
inventory control with varying demand, and project scheduling with random task
duration, etc. While the effect of uncertain parameters on the optimal
solution and objective value can be studied through the well-known sensitivity
analysis, or what-if analysis, in a deterministic optimisation approach, such
post-optimality analysis does not prescribe solutions under uncertainty a
priori. This subsection provides an overview of a suite of optimisation models
and methods that seek to obtain optimal or near-optimal solutions for the
class of decision problems where some parameters are stochastic with known
probability distribution373737Decision problems with incomplete information
about random parameter’s probability distribution go beyond the scope of this
review..

Originated in the seminal work of Dantzig (1955), stochastic programming is
one of the earliest and most prominent approaches to deal with optimisation
problems with stochastic parameters. The basic stochastic programming model
has a two-stage framework, called two-stage stochastic programming with
recourse (2S-SPR; Birge and Louveaux, 2011). In the first stage, the here-and-
now decision is made. Then in the second stage, the recourse decision is
prescribed for each scenario of stochastic parameters after their realisation.
The objective function minimises the total cost as the summation of the first-
stage cost and the expected second-stage cost given the probability
distribution of the stochastic parameters. It is often insightful to compute
the value of stochastic solution (VSS; Birge, 1982) as the difference between
the optimal objective function value of the deterministic counterpart (by
substituting the stochastic parameters with their point estimates) and that of
the stochastic programming model. We refer to Birge and Louveaux (2011) and
Shapiro et al. (2021) for a systematic and updated treatment on the modelling
and theory of stochastic programming, and to Wallace and Ziemba (2005) for a
collection of applications of stochastic programming. Recent applications
include disaster relief management (Grass and Fischer, 2016), transit network
design (Zhao et al., 2017), portfolio selection (Masmoudi and Abdelaziz,
2018), treatment plant placement in drinking water systems (Schwetschenau et
al., 2019), process systems (Li and Grossmann, 2021), multi-product aggregate
planning (Gómez-Rocha and Hernández-Gress, 2022), and resource allocation for
infrastructure planning (Zhang and Alipour, 2022), among others.

A stochastic programming model may also include a constraint that is satisfied
with a probability. This model is is known as the chance constrained
programming model introduced by Charnes et al. (1958). The probabilistic
constraint can often be transformed into a deterministic constraint given the