Journal of the Operational Research Society

ISSN: (Print) (Online) Journal homepage: https://www.tandfonline.com/loi/tjor20

Operational Research: methods and applications

Fotios Petropoulos, Gilbert Laporte, Emel Aktas, Sibel A. Alumur, Claudia Archetti, Hayriye Ayhan, Maria Battarra, Julia A. Bennell, Jean-Marie Bourjolly, John E. Boylan, Michèle Breton, David Canca, Laurent Charlin, Bo Chen, Cihan Tugrul Cicek, Louis Anthony Cox Jr, Christine S.M. Currie, Erik Demeulemeester, Li Ding, Stephen M. Disney, Matthias Ehrgott, Martin J. Eppler, Güneş Erdoğan, Bernard Fortz, L. Alberto Franco, Jens Frische, Salvatore Greco, Amanda J. Gregory, Raimo P. Hämäläinen, Willy Herroelen, Mike Hewitt, Jan Holmström, John N. Hooker, Tuğçe Işık, Jill Johnes, Bahar Y. Kara, Özlem Karsu, Katherine Kent, Charlotte Köhler, Martin Kunc, Yong-Hong Kuo, Adam N. Letchford, Janny Leung, Dong Li, Haitao Li, Judit Lienert, Ivana Ljubić, Andrea Lodi, Sebastián Lozano, Virginie Lurkin, Silvano Martello, Ian G. McHale, Gerald Midgley, John D.W. Morecroft, Akshay Mutha, Ceyda Oğuz, Sanja Petrovic, Ulrich Pferschy, Harilaos N. Psaraftis, Sam Rose, Lauri Saarinen, Said Salhi, Jing-Sheng Song, Dimitrios Sotiros, Kathryn E. Stecke, Arne K. Strauss, İstenç Tarhan, Clemens Thielen, Paolo Toth, Tom Van Woensel, Greet Vanden Berghe, Christos Vasilakis, Vikrant Vaze, Daniele Vigo, Kai Virtanen, Xun Wang, Rafał Weron, Leroy White, Mike Yearworth, E. Alper Yıldırım, Georges Zaccour & Xuying Zhao

**To cite this article:** Fotios Petropoulos, Gilbert Laporte, Emel Aktas, Sibel A. Alumur, Claudia Archetti, Hayriye Ayhan, Maria Battarra, Julia A. Bennell, Jean-Marie Bourjolly, John E. Boylan, Michèle Breton, David Canca, Laurent Charlin, Bo Chen, Cihan Tugrul Cicek, Louis Anthony Cox Jr, Christine S.M. Currie, Erik Demeulemeester, Li Ding, Stephen M. Disney, Matthias Ehrgott, Martin J. Eppler, Güneş Erdoğan, Bernard Fortz, L. Alberto Franco, Jens Frische, Salvatore Greco, Amanda J. Gregory, Raimo P. Hämäläinen, Willy Herroelen, Mike Hewitt, Jan Holmström, John N. Hooker, Tuğçe Işık, Jill Johnes, Bahar Y. Kara, Özlem Karsu, Katherine Kent, Charlotte Köhler, Martin Kunc, Yong-Hong Kuo, Adam N. Letchford, Janny Leung, Dong Li, Haitao Li, Judit Lienert, Ivana Ljubić, Andrea Lodi, Sebastián Lozano, Virginie Lurkin, Silvano Martello, Ian G. McHale, Gerald Midgley, John D.W. Morecroft, Akshay Mutha, Ceyda Oğuz, Sanja Petrovic, Ulrich Pferschy, Harilaos N. Psaraftis, Sam Rose, Lauri Saarinen, Said Salhi, Jing-Sheng Song, Dimitrios Sotiros, Kathryn E. Stecke, Arne K. Strauss, İstenç Tarhan, Clemens Thielen, Paolo Toth, Tom Van Woensel, Greet Vanden Berghe, Christos Vasilakis, Vikrant Vaze, Daniele Vigo, Kai Virtanen, Xun Wang, Rafał Weron, Leroy White, Mike Yearworth, E. Alper Yıldırım, Georges Zaccour & Xuying Zhao (27 Dec 2023): Operational Research: methods and applications, Journal of the Operational Research Society, DOI: 10.1080/01605682.2023.2253852

**To link to this article:** https://doi.org/10.1080/01605682.2023.2253852

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Submit your article to this journal

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

View related articles

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

CrossMark

View Crossmark data

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

REVIEW ARTICLE

<figure><img src="image_11.png" /></figure>

OPEN ACCESS

<figure><img src="image_12.png" /></figure>

Check for updates

Operational Research: methods and applications

Fotios Petropoulos<sup>a,b</sup>, Gilbert Laporte<sup>a,c,d</sup>, Emel Aktas<sup>e</sup> <sup>ID</sup>, Sibel A. Alumur<sup>f</sup>, Claudia Archetti<sup>g</sup>, Hayriye Ayhan<sup>h</sup>, Maria Battarra<sup>a</sup>, Julia A. Bennell<sup>i</sup>, Jean-Marie Bourjolly<sup>j</sup>, John E. Boylan<sup>k</sup> <sup>ID</sup>, Michèle Breton<sup>c</sup>, David Canca<sup>l</sup> <sup>ID</sup>, Laurent Charlin<sup>c,m</sup>, Bo Chen<sup>n</sup> <sup>ID</sup>, Cihan Tugrul Cicek<sup>o</sup>, Louis Anthony Cox Jr<sup>p,q</sup>, Christine S.M. Currier<sup>r</sup>, Erik Demeulemeester<sup>s</sup>, Li Ding<sup>t</sup>, Stephen M. Disney<sup>u</sup>, Matthias Ehrgott<sup>v</sup>, Martin J. Eppler<sup>w</sup>, Güneş Erdoğan<sup>a</sup>, Bernard Fortz<sup>x,y,z</sup>, L. Alberto Franco<sup>aa,ab</sup>, Jens Frische<sup>ac</sup>, Salvatore Greco<sup>ad,ae</sup>, Amanda J. Gregory<sup>af</sup>, Raimo P. Hämäläinen<sup>ag</sup>, Willy Herroelen<sup>s</sup>, Mike Hewitt<sup>ah</sup>, Jan Holmström<sup>ai</sup>, John N. Hooker<sup>aj</sup>, Tuğçe Işık<sup>ak</sup>, Jill Johnes<sup>al</sup>, Bahar Y. Kara<sup>am</sup>, Özlem Karsu<sup>am</sup> <sup>ID</sup>, Katherine Kent<sup>an</sup>, Charlotte Köhler<sup>ao</sup>, Martin Kunc<sup>ap</sup>, Yong-Hong Kuo<sup>aq,ar</sup>, Adam N. Letchford<sup>v</sup> <sup>ID</sup>, Janny Leung<sup>as</sup>, Dong Li<sup>aa</sup>, Haitao Li<sup>at</sup> <sup>ID</sup>, Judit Lienert<sup>au</sup>, Ivana Ljubić<sup>av</sup>, Andrea Lodi<sup>aw</sup>, Sebastián Lozano<sup>ax</sup> <sup>ID</sup>, Virginie Lurkin<sup>ay</sup>, Silvano Martello<sup>az</sup>, Ian G. McHale<sup>ba</sup>, Gerald Midgley<sup>af,bb,bc,bd,be*</sup>, John D.W. Morecroft<sup>bf</sup>, Akshay Mutha<sup>b9</sup>, Ceyda Oğuz<sup>bh</sup>, Sanja Petrovic<sup>bi</sup>, Ulrich Pferschy<sup>bj</sup>, Harilaos N. Psaraftis<sup>bk</sup>, Sam Rose<sup>bl</sup>, Lauri Saarinen<sup>ai</sup>, Said Salhi<sup>bm</sup>, Jing-Sheng Song<sup>bn</sup>, Dimitrios Sotiros<sup>bo</sup>, Kathryn E. Stecke<sup>bp</sup>, Arne K. Strauss<sup>ac</sup>, İstenç Tarhan<sup>bh</sup>, Clemens Thielen<sup>bq</sup>, Paolo Toth<sup>az</sup>, Tom Van Woensel<sup>by</sup>, Greet Vanden Berghe<sup>br</sup>, Christos Vasilakis<sup>a</sup>, Vikrant Vaze<sup>bs</sup>, Daniele Vigo<sup>bt</sup>, Kai Virtanen<sup>bu,bv</sup>, Xun Wang<sup>bw</sup>, Rafał Weron<sup>bo</sup>, Leroy White<sup>bx</sup>, Mike Yearworth<sup>u</sup>, E. Alper Yıldırım<sup>bz</sup>, Georges Zaccour<sup>c</sup> and Xuying Zhao<sup>ca</sup>

<sup>a</sup>School of Management, University of Bath, Bath, UK; <sup>b</sup>Makridakis Open Forecasting Center, University of Nicosia, Nicosia, Cyprus; <sup>c</sup>Department of Decision Sciences, HEC Montréal, Montreal, Canada; <sup>d</sup>Molde University College, Molde, Norway; <sup>e</sup>Cranfield School of Management, Cranfield University, Cranfield, UK; <sup>f</sup>Department of Management Sciences, University of Waterloo, Waterloo, Canada; <sup>g</sup>Department of Information Systems, Decision Sciences and Statistics, ESSEC Business School in Paris, Cergy, France; <sup>h</sup>H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, USA; <sup>i</sup>Centre for Decision Research, Leeds University Business School, University of Leeds, Leeds, United Kingdom; <sup>j</sup>Université du Québec à Montréal, Montreal, Canada; <sup>k</sup>Centre for Marketing Analytics and Forecasting, Lancaster University Management School, Lancaster University, Lancaster, UK; <sup>l</sup>Department of Industrial Engineering and Management Science I, School of Engineering, University of Seville, Seville, Spain; <sup>m</sup>Mila-Quebec AI Institute, Montreal, Canada; <sup>n</sup>Warwick Business School, University of Warwick, Coventry, UK; <sup>o</sup>Department of Industrial Engineering, Atilim University, Ankara, Turkey; <sup>p</sup>Department of Business Analytics, University of Colorado School of Business, Denver, USA; <sup>q</sup>Cox Associates, Denver, USA; <sup>r</sup>Mathematical Sciences Department, University of Southampton, Southampton, UK; <sup>s</sup>Faculty of Economics and Business, Research Center for Operations Management, KU Leuven, Leuven, Belgium; <sup>t</sup>Durham University Business School, Durham University, Durham, UK; <sup>u</sup>Center for Simulation, Analytics, and Modelling, University of Exeter Business School, Exeter, UK; <sup>v</sup>Department of Management Science, Lancaster University Management School, Lancaster University, Lancaster, UK; <sup>w</sup>University of St. Gallen, St. Gallen, Switzerland; <sup>x</sup>HEC – Management School of the University of Liège, Liège, Belgium; <sup>y</sup>Département d'Informatique, Université Libre de Bruxelles, Brussels, Belgium; <sup>z</sup>Inria Lille-Nord Europe, Villeneuve d'Ascq, France; <sup>aa</sup>School of Business and Economics, Loughborough University, Loughborough, UK; <sup>ab</sup>Universidad del Pacifico, Lima, Perú; <sup>ac</sup>WHU – Otto Beisheim School of Management, Vallendar, Germany; <sup>ad</sup>Department of Economics and Business, University of Catania, Catania, Italy; <sup>ae</sup>Portsmouth Business School, Centre of Operations Research and Logistics (CORL), University of Portsmouth, Portsmouth, UK; <sup>af</sup>Faculty of Business, Law and Politics, Centre for Systems Studies, University of Hull, Hull, UK; <sup>ag</sup>Systems Analysis Laboratory, Aalto University, Finland; <sup>ah</sup>Department of Information Systems and Supply Chain Management, Loyola University, Chicago, USA; <sup>ai</sup>Department of Industrial Engineering and Management, Aalto University, Aalto, Finland; <sup>aj</sup>Tepper School of Business, Carnegie Mellon University, Pittsburgh, USA; <sup>ak</sup>Department of Industrial Engineering, Clemson University, Clemson, USA; <sup>al</sup>Department of Accounting, Finance and Economics, University of Huddersfield, Huddersfield, UK; <sup>am</sup>Department of Industrial Engineering, Bilkent University, Ankara, Turkey; <sup>an</sup>Office for National Statistics, Newport, UK; <sup>ao</sup>Europa Universität Viadrina, Frankfurt (Oder), Germany; <sup>ap</sup>Southampton Business School, University of Southampton, Southampton, UK; <sup>aq</sup>Department of Industrial and Manufacturing Systems Engineering, The University of Hong Kong, Pokfulam, Hong Kong, China; <sup>at</sup>HKU Musketeers Foundation Institute of Data Science, The University of Hong Kong, Pokfulam, Hong Kong, China; <sup>as</sup>State Key Laboratory of Internet of Things for Smart City, Choi Kai Yau College, University of Macau, Taipa, Macau, China; <sup>at</sup>Supply Chain & Analytics Department, University of Missouri-St. Louis, St. Louis, USA; <sup>au</sup>Eawag: Swiss Federal Institute of Aquatic Science and Technology, Dübendorf, Switzerland; <sup>av</sup>ESSEC Business School, Cergy-Pontoise, France; <sup>aw</sup>Jacobs Technion-Cornell Institute, Cornell Tech and Technion – IIT, New York, USA; <sup>ax</sup>Department of Industrial Management, Escuela Superior de Ingenieros, University of Seville, Sevilla, Spain; <sup>ay</sup>HEC Lausanne Faculty of Business and Economics, University of Lausanne, Lausanne, Switzerland; <sup>az</sup>DEI “Guglielmo Marconi”, Alma Mater Studiorum Università di Bologna, Bologna, Italy; <sup>ba</sup>Centre for Sports Business, University of Liverpool Management School, UK; <sup>bb</sup>Birmingham Leadership Institute, University of Birmingham, Birmingham, UK; <sup>bc</sup>Department of Informatics, Faculty of Technology, Linnaeus University, Växjö, Sweden; <sup>bd</sup>School of Innovation, Design and Engineering, Mälardalen University, Eskilstuna, Sweden; <sup>be</sup>School of Agriculture and Food Sciences, University of Queensland, Brisbane, Queensland, Australia; <sup>bf</sup>London Business School, London, UK; <sup>bg</sup>Grossman School of Business, University of Vermont, Burlington, Vermont, USA; <sup>bh</sup>Department of Industrial Engineering, Koç University, İstanbul, Turkey; <sup>bl</sup>Nottingham University Business School, University of Nottingham, Nottingham, UK; <sup>bj</sup>Department of

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>

Operations and Information Systems, University of Graz, Graz, Austria; ᵇDepartment of Technology, Management and Economics, Technical University of Denmark, Lyngby, Denmark; ᵇDepartment for Transport, London, UK; ᵇCentre for Logistics and Heuristic optimisation, Kent Business School, University of Kent, Kent, UK; ᵇThe Fuqua School of Business, Duke University, Durham, North Carolina, USA; ᵇDepartment of Operations Research and Business Intelligence, Wrocław University of Science and Technology, Wrocław, Poland; ᵇNaveen Jindal School of Management, The University of Texas at Dallas, Richardson, Texas, USA; ᵇTUM Campus Straubing for Biotechnology and Sustainability, Technical University of Munich, Straubing, Germany; ᵇSchool of Industrial Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands; ᵇDepartment of Computer Science, KU Leuven, Gent, Belgium; ᵇThayer School of Engineering, Dartmouth College, Hanover, USA; ᵇDepartment of Electrical, Electronic, and Information Engineering "G. Marconi" and CIRI-ICT, University of Bologna, Bologna, Italy; ᵇDepartment of Mathematics and Systems Analysis, Aalto University, Helsinki, Finland; ᵇDepartment of Military Technology, National Defence University, Helsinki, Finland; ᵇLogistics Systems Dynamics Group, Cardiff Business School, Cardiff University, UK; ᵇManagement Department, University of Exeter Business School, Exeter, UK; ᵇSchool of Mathematics, The University of Edinburgh, Edinburgh, UK; ᵈMendoza College of Business, University of Notre Dame, Notre Dame, Indiana, USA

## ABSTRACT

Throughout its history, Operational Research has evolved to include methods, models and algorithms that have been applied to a wide range of contexts. This encyclopedic article consists of two main sections: methods and applications. The first summarises the up-to-date knowledge and provides an overview of the state-of-the-art methods and key developments in the various subdomains of the field. The second offers a wide-ranging list of areas where Operational Research has been applied. The article is meant to be read in a nonlinear fashion and used as a point of reference by a diverse pool of readers: academics, researchers, students, and practitioners. The entries within the methods and applications sections are presented in alphabetical order.

The authors dedicate this paper to the 2023 Turkey/Syria earthquake victims. We sincerely hope that advances in OR will play a role towards minimising the pain and suffering caused by this and future catastrophes.

Operations research is neither a method nor a technique; it is or is becoming a science and as such is defined by a combination of the phenomena it studies.

Ackoff (1956)

# 1. Introduction¹

The year 2024 marks the 75th anniversary of the *Journal of the Operational Research Society*, formerly known as *Operational Research Quarterly*. It is the oldest Operational Research (OR) journal worldwide. On this occasion, my colleague Fotios Petropoulos from University of Bath proposed to the editors of the journal to edit an encyclopedic article on the state of the art in OR. Together, we identified the main methodological and application areas to be covered, based on topics included in the major OR journals and conferences. We also identified potential authors who responded enthusiastically and whom we thank wholeheartedly for their contributions.

Modern OR originated in the United Kingdom during World War II as a need to support the operations of early radar-detecting systems and was later applied to other operations (McCloskey, 1987). However, one could argue that it precedes this period in history since it is partly rooted in several mathematical fields such as probability theory and statistics, calculus, and linear algebra, developed much earlier. For example, the Fourier-Motzkin elimination method (Fourier, 1826a, 1826b) constitutes the main basis of linear programming. Queueing

theory, which plays a central role in telecommunications and computing, already existed as a distinct field of study since the early 20th century (Erlang, 1909), and other concepts, such as the economic order quantity (Harris, 1913) were developed more than one century ago. Interestingly, while many recent advances in OR are rooted in theoretical or algorithmic concepts, we are now witnessing a return to the practical roots of OR through the development of new disciplines such as business analytics.

After the war ended, several industrial applications of OR arose, particularly in the manufacturing and mining sectors which were then going through a renaissance. The transportation sector is without doubt the field that has most benefited from OR, mostly since the 1960s. The aviation, rail, and e-commerce industries could simply not operate at their current scale without the support of massive data analysis and sophisticated optimisation techniques. The application of OR to maritime transportation is more recent, but it is fast gaining in importance. Other areas that are less visible, such as telecommunications, also deeply depend on OR. The success of OR in these fields is partly explained by their network structures which make them amenable to systematic analysis and treatment through mathematical optimisation techniques. In the same vein, OR also plays a major role in various branches of logistics and project management, such as facility location, forecasting, inventory planning, scheduling, and supply chain management.

The public sector and service industries also benefit greatly from OR. Healthcare is the first area that

## ARTICLE HISTORY

Received 24 March 2023

Accepted 26 August 2023

## KEYWORDS

Review; encyclopedia;

theory; practice; principles;

optimisation; programming;

systems; simulation;

decision making; models
<figure><img src="image_16.png" /></figure>

comes to mind because of its very large scale and complexity. Decision making in healthcare is more decentralised than in transportation and manufacturing, for example, and the human issues involved in this sector add a layer of complexity. OR methodologies have also been applied to diverse areas such as education, sports management, natural resources, environment and sustainability, political districting, safety and security, energy, finance and insurance, revenue management, auctions and bidding, and disaster relief, most of which are covered in this article.

Among OR methodologies, mathematical programming occupies a central place. The simplex method for linear programming, conceived by Dantzig in 1947 but apparently first published later (Dantzig, 1951), is arguably the single most significant development in this area. Over time, linear programming has branched out into several fields such as nonlinear programming, mixed integer programming, network optimisation, combinatorial optimisation, and stochastic programming. The techniques most frequently employed for the exact solution of mathematical programs are based on branch-and-bound, branch-and-cut, branch-and-price (column generation), and dynamic programming. Game theory and data envelopment analysis are firmly rooted in mathematical programming. Control theory is also part of continuous mathematical optimisation and relies heavily on differential equations.

Complexity theory is fundamental in optimisation. Most problems arising in combinatorial optimisation are NP-hard and typically require the application of heuristics for their solution. Much progress has been made in the past 40 years or so in the development of metaheuristics based on local search, genetic search, and various hybridisation schemes. Many problems in fields such as vehicle routing, location analysis, cutting and packing, set covering, and set partitioning can now be solved to near optimality for realistic sizes by means of modern heuristics. A recent trend is the use of open-source software which not only helps disseminate research results, but also contributes to ensuring their accuracy, reproducibility and adoption.

Several modelling paradigms such as systems thinking and systems dynamics approach problems from a high-level perspective, examining the inter-relationships between multiple elements. Complex systems can often be analysed through simulation, which is also commonly used to assess the performance of heuristics. Decision analysis provides a useful framework for structuring and solving complex problems involving soft and hard criteria, behavioural OR, stochasticity, and dynamism. Recently, issues related to ethics and fairness have come to play an increasing role in decision making.

Because the various topics of this review paper are listed in alphabetical order, the subsection on “Artificial intelligence, machine learning and data science” comes first, but this topic constitutes one of the latest developments in the field. It holds great potential for the future and is likely to reshape parts of the OR discipline. Already, machine learning-based heuristics are competitive for the solution of some hard problems.

This paper begins with a quote from Russell L. Ackoff who has been a pioneer of OR. In 1979, he published in this journal two articles (Ackoff, 1979a, 1979b) that presented a rather pessimistic view of our discipline. The author complained about the lack of communications between academics and practitioners, and about the fact that some OR curricula in universities did not sufficiently prepare students for practice, which is still true to some extent. One of his two articles is entitled “The Future of Operational Research is Past”, which may be perceived as an overreaction to this diagnosis. In my view, the present article provides clear evidence to the contrary. Soon after the publication of the two Ackoff papers, we witnessed the development of micro-computing, the Internet and the World Wide Web. It has become much easier for researchers in our community to access information, software and computing facilities, and for practitioners to access and use our research results. We are now fortunate to have access to sophisticated open-source software, data bases, bibliographic sources, editing and visualisation tools, and communication facilities. Our field is richer than it has ever been, both in terms of theory and applications. It is constantly evolving in interaction with other disciplines, and it is clearly alive and well and has a promising future.

## 2. Methods

### 2.1. Artificial intelligence, machine learning, and data science²

Machine learning (ML) comprises techniques for modelling predictive tasks, i.e., tasks that involve the prediction of an unknown quantity from other observed quantities. Ideas of learning in an artificial system and the term machine learning were first discussed in the 1950s (Samuel, 1959) and their development and popularity have seen enormous growth over the last two decades in part due to the availability of large-scale datasets and increased computational resources to model them.

Mitchell (1997) provides this concrete definition of machine learning “A computer program is said to learn from experience E with respect to some class of tasks T, and performance measure P, if its performance at tasks in T, as measured by P, improves with experience
<figure><img src="image_17.png" /></figure>

E". The program is a model or a function and its experience E is the type of data it has access to. There are three types of experiences supervised, unsupervised, and reinforcement learning. The performance measure (P) allows for model evaluation and comparison including model selection.

Supervised learning is an experience where a model aims at predicting one or more unobserved target (dependent) variables given observed ibackpr-out (independent) variables. In other words, a super-vised model is a function that map inputs to outputs. The process of solving a supervised problem involves first learning a model, that is adjusting its parameters using a training dataset with both input and target variables. The training set is drawn IID (independ-ently and identically distributed) from an underlying distribution over inputs and targets. Once trained, the model can provide target predictions for new unseen samples from the same distribution. The most common tasks in supervised learning are regression (real dependent variable) and classification (categorical dependent variable). Evaluating a super-vised system is usually performed using held-out data referred to as the test data while held-out validation data is used for model development and selection using procedures such as k-fold cross-validation.

Supervised models can be dichotomised into linear and nonlinear models. Linear models perform a linear mapping from inputs to outputs (e.g., linear regression). Machine learning mostly investigates nonlinear super-vised models including deep neural network (DNN) models (Goodfellow et al., 2016). DNNs are composed of a succession of parametrised nonlinear transforma-tions called layers and each layer contains a set of transformations called neurons. Layers successively transform an input datum into a target. The parameters of the layers are adjusted to iteratively obtain better predictions using a procedure called backpropagation, a form of gradient descent (Goodfellow et al., 2016, §6.5). DNNs are state-of-the-art methods for many large-scale non-structured datasets across domains (see also §3.19). DNNs can be adapted to different sizes of inputs and targets as well as variable types. They can also be speci-fised for specific types of data. Recurrent neural net-works (RNNs) are auto-regressive models for sequential data (Rumelhart et al., 1986). The sequential data are tokenised and an RNN transforms each token sequen-tially along with a transformation of the previous tokens. Convolutional neural networks (CNNs) are spe-cialised networks for modelling data that is arranged on a grid (e.g., an image Lecun, 1989). Their layers contain a convolution operation between an input and a parameterised filter followed by a nonlinear trans-formation, and a pooling operation. Each layer proc-esses data locally and so requires fewer parameters compared to vanilla DNNs. As a result, CNNs can

model higher-dimensional data. Graphical neural net-
works (GNNs) are specialised architectures for model-
ling graph data (e.g., a social network; Scarselli et al.,
2009). In GNNs, the data are transformed by following
the topology of the graph. Last, attention layers dynam-
ically combine their inputs (tokens) based on their val-
ues. Transformer models use successions of attention
and feed-forward layers to model sequential input and
output data (Vaswani et al., 2017). Transformers are
more efficient to train than RNNs and can be trained
on internet-scale data given enormous computational
power. The availability of such broad datasets especially
in the text and image domains has given rise to a class
of very-large-scale models (also referred to as founda-
tion models) that display an ability to adapt to and
obtain high performance across a diversity of down-
stream supervised tasks (Bommasani et al., 2021)

Last, attention is a mechanism that considers data
to be unordered and uses transformations dynamic-
ally. Transformers are models based on attention.
They provide more efficient training than RNNs for
very large-scale datasets (Vaswani et al., 2017).

Neural networks currently outperform other meth-
ods when learning from unstructured data (e.g.,
images and text). For tabular data, data that is natur-
ally encoded in a table and that has heterogeneous
features (Grinsztajn et al., 2022), best-performing
methods use ideas first proposed in tree-based classi-
fiers, bagging, and boosting. They include random
forests (Breiman, 2001), XGBoost (Chen & Guestrin,
2016) which both scale to large-scale datasets as well
as kernel methods including support vector machines
(SVMs see, e.g., Schlkopf et al., 2018) and probabilis-
tic Gaussian Processes (GPs see, e.g., Rasmussen &
Williams, 2005). These methods are used across
regression and classification tasks.

In unsupervised learning, the second type of
experience, the data consist of independent variables
 fea tures or covariates) alone. The aim of unsupe-
 r vised learning is to model the structure of the data
to better understand their properties. As a result,
evaluating an unsupervised model is often task and
application-dependant (Murphy, 2022, §1.3.4). The
prototypical unsupervised-learning task is clustering.
It involves learning a function that groups similar
data together according to a similarity measure and
desiderata often expressed as an objective function.
Several standard algorithms divided into hierarchical
and non-hierarchical methods exist. The former
uses the similarity between all pairs of data and
finds a hierarchy of clustering solutions with a dif-
ferent number of clusters using either a bottom-up
or top-down approach. Agglomerative clustering is a
standard hierarchical approach. Non-hierarchical
methods tend to be more computationally efficient
in terms of dataset size. For example, K-means
<figure><img src="image_18.png" /></figure>

clustering is a well-known non-hierarchical method that finds a single solution using *K* clusters (MacQueen, 1967). Other unsupervised learning tasks include dimensionality reduction for example for visualisation or to prepare data for further analysis. Density modelling is another unsupervised task where a probabilistic model learns to assign a probability to each datum (Murphy, 2022, §1.3). Probabilistic models can be used to learn the hidden structure in large quantities of data (e.g., Hoffman et al., 2013). Further, probabilistic models are also used to generate high-dimensional data (e.g., images of human faces or English text) with high fidelity (Karras et al., 2021) and often referred in this context as generative models. Large Language Models are examples of such generative models (Bommasani et al., 2021).

Reinforcement learning (RL) is the third type of experience. RL models collect their own data by executing actions in their environment to maximise their reward. RL is a sequential decision-making task and is formalised using Markov decision processes (MDPs) (Sutton & Barto, 2018, §3.8). An MDP encodes a set of states, available actions, distribution over next states given current states and action, a reward function, and a discount factor. Partially observable MDPs (or POMDPs) extend the formalism to environments where the exact current state is unknown (Kaelbling et al., 1998). In RL, an agent's objective is to learn a policy, a distribution over actions for each state in an environment. Tasks are defined by rewards attached to different states. Exact and approximate methods exist for solving RL problems. Whereas exact solutions are appropriate for smaller tabular problems only, deep neural networks are widely used for solving larger-scale problems that require approximate solutions yielding a set of techniques known as deep reinforcement learning (Mnih et al., 2015). An RL agent can also learn to imitate an expert either by learning a mapping from states/observations to actions as in supervised learning (a technique known as imitation learning; for a survey, see Hussein et al., 2017) or by trying to learn the expert's reward function (inverse reinforcement learning Russell, 1998).

In addition to learning models for solving prediction tasks using one of the three experiences above, machine learning also studies methods for enabling the reuse of information learned from one or multiple datasets and environments to other similar ones. Representation learning studies how to learn such reusable information and it can use both supervised and unsupervised experiences (Murphy, 2023, §32). When using a deep learning model, a representation is obtained after one or more layer transformations of the data. Representation learning

is used in a variety of situations including for transfer learning tasks, where a trained model is reused to solve a different supervised task (for a survey, see Zhuang et al., 2021).

In the last decade, machine learning models have achieved high performance on a variety of tasks including perceptual ones (e.g., recognising objects in images and words from speech) as well as natural language processing ones thereby becoming a core component of artificial intelligence (AI) methods. The goal of AI methods is to develop intelligent systems. Some of these advances shine a bright light on the ethical aspects of machine learning techniques and are active areas of study (see, e.g., Dignum, 2019; Barocas et al., 2019). Another area of active study is explainability (Phillips et al., 2021). Some of the most effective ML tools make predictions and recommendations that are hard to explain to users (for example when neural networks are employed). Clearly, lack of explainability slows down ML use in those contexts where decisions made due to those predictions and recommendations are life changing and involve a human in the loop, healthcare (applying a treatment), finance (refusing a mortgage), or justice (granting parole) to mention a few. So, explainability is currently one of the most crucial challenges for ML and AI and, at the same time, a tremendous opportunity for their wider applicability.

Further, advances in machine learning alongside statistics, data management, and data processing, as well as the wider availability of datasets from a variety of domains have led to the popularisation and development of data science (DS), a discipline whose goal is to extract insights and knowledge from these data. DS uses statistics and machine-learning techniques for inference and prediction, but it also aims at enabling and systematising the analysis of large quantities of data. As such, it includes components of data management, visualisation, as well as the design of (efficient) data processing algorithms (Grus, 2019).

### 2.1.1. Resources

Murphy (2022) provides a thorough introduction to the field following a probabilistic approach and its sequel (Murphy, 2023) introduces advanced topics. Goodfellow et al. (2016) provide a self-contained introduction to the field of deep learning (the field evolves rapidly and more advanced topics are covered through recent papers and in Murphy, 2023). Open-source software packages in Python and other languages are essential. They include data-wrangling libraries such as pandas (McKinney, 2010) and plotting ones such as matplotlib (Hunter, 2007). The library scikit-learn (Pedregosa et al., 2011) in Python offers an extensive API that includes data
<figure><img src="image_19.png" /></figure>

processing, a toolbox of standard supervised and unsupervised models, and evaluation routines. For deep learning, PyToch (Paszke et al., 2019) and TensorFlow (Abadi et al., 2015) are the standard.

### 2.1.2. Learning for combinatorial optimisation

The impressive success of machine learning in the last decade made it natural to explore its use in many scientific disciplines, such as drug discovery and material sciences. Combinatorial optimisation (CO; §2.4) is no exception to this trend and we have witnessed an intense exploration (or, better, revival) of the use of machine learning for CO. Two lines of work have strongly emerged. On the one side, ML has been used to learn crucial decisions within CO algorithms and solvers. This includes imitating an algorithmic expert that is computationally expensive like in the case of strong branching for branch and bound, the single application that has attracted the largest amount of interest (Lodi & Zarpellon, 2017; Gasse et al., 2019). The interested reader is referred to two recent surveys (Bengio et al., 2021; Cappart et al., 2021), the latter highlighting the relevance of GNNs for effective CO representation. On the other side, ML has been used end to end, i.e., for solving CO problems directly or leveraging ML to devise hybrid methods for CO. The area is surveyed in Kotary et al. (2021).

### 2.2. Behavioural OR³

Behavioural OR (BOR) is concerned with the study of human behaviour in OR-supported settings. Specifically, BOR examines how the behaviour of individuals affects, or is affected by, an OR-supported intervention⁴. The individuals of interest are those who, acting in isolation or as part of a team, design, implement and engage with OR in practice. These individuals include OR practitioners playing specific intervention roles (e.g., modellers, facilitators, consultants), and other individuals with varying interests and stakes in the intervention (e.g., users, clients, domain experts, sponsors).

A concern with the behavioural aspects of the OR profession can be traced back to past debates in the 1960s, 1970s and 1980s (Churchman, 1970; Dutton & Walton, 1964; Jackson et al., 1989). Although these debates dwindled down in subsequent years, the emergence of BOR as a field of study represents a return to these earlier concerns (Franco & Hämäläinen, 2016; Hämäläinen et al., 2013). What motivates this resurgence is the recognition that the successful deployment of OR in practice relies heavily on our understanding of human behaviour. For example, overconfidence, competing interests, and the willingness to expend effort in

searching, sharing, and processing information are three behavioural issues that can negatively affect the success of OR activities. Attention to behavioural issues has been central in disciplines such as economics, psychology and sociology for decades, and BOR studies draw heavily from these reference disciplines (Franco et al., 2021).

It is important to distinguish between the specific focus of BOR and the broader focus of behavioural modelling. The creation of models that capture human behaviour has a long tradition within OR, but it is not necessarily concerned with the study of human behaviour in OR-supported settings. For example, in the last 20 years operational researchers have produced an increasing number of robust analytical models that describe behaviour in, and predict its impact on, operations management settings (Cui & Wu, 2018; Donohue et al., 2020; Loch & Wu, 2007). Operational researchers also have produced simulation models that capture human behaviour within a system with different levels of complexity. For example, systems dynamics models incorporate high-level variables representing average behaviour (Morecroft, 2015; Sterman, 2000, §2.22), and discrete event simulation models capture human processes controlled by simple behavioural rules (Brailsford & Schmidt, 2003; Robinson, 2014, §2.19). More complex agent-based simulation models represent behaviour as emergent from the interactions of agents with particular behavioural attributes (Sonnessa et al., 2017; Utomo et al., 2018, §2.19). Overall, behavioural modelling within the OR field is concerned with examining human behaviour in a system of interest in order to improve that system⁵. In contrast, BOR takes an OR-supported intervention as the core system of interest where human behaviour is examined. The ultimate goal of BOR is to generate an improved understanding of the behavioural dimension of OR practice, and use this understanding to design and implement better OR-supported interventions.

Another important distinction worth stating is that between BOR and Soft OR. At first glance, this distinction may seem unnecessary as BOR is a field of study within OR, while Soft OR refers to a specific family of problem structuring approaches (§2.20). Soft OR approaches have been developed to help groups reach agreements on problem structure and, often, appropriate responses to a problem of concern (Franco & Rouwette, 2022; Rosenhead & Mingers, 2001). However, while Soft OR intervention design and implementation typically require the consideration of behavioural issues, this is not the same as choosing human behaviour in a Soft OR intervention context as the unit of analysis. Of course, a study with such a focus would certainly
<figure><img src="image_20.png" /></figure>

fall within the BOR remit (e.g., Tavella et al., 2021). But note that BOR is also concerned with the study of human behaviour in other OR-supported settings, such as those involving the use of 'hard' and 'mixed-method' OR approaches.

Studies of behaviour in OR-supported settings assume implicitly or explicitly that human behaviour is either influenced by cognitive and external factors, or is in itself an influencing factor (Franco et al., 2021). In the first case, observed individual and collective action is taken to be guided by cognitive structures (e.g., personality traits, cognitive styles) manifested during OR-related activity – behaviour is *influenced*. In contrast, the second case assumes that individuals and collectives are responsible for determining how OR-related activity will unfold – behaviour is *influencing*. This raises the practical possibility that the same OR methodology, technique, or model could be used in distinctive ways by various individuals or groups according to their cognitive orientations, goals and interests (Franco, 2013). Whilst behaviour in practice is likely to lie somewhere between the influenced and influencing assumptions, BOR studies tend to foreground one of the extremes as the focus, while backgrounding the other.

BOR studies can adopt three different research methodologies to examine behaviour: variance, process, and modelling. A *variance methodology* uses variables that represent the important aspects or attributes of the OR-supported activity being examined. Variance explanations of behavioural-related phenomena take the form of causal statements captured in a theoretically-informed research model that incorporates these variables (e.g., A causes B, which causes C). The research model is then tested with data generated by the activity, and the research findings are assessed in terms of their generality (Poole, 2004). Adopting a variance research methodology typically requires the implementation of experimental, quasi-experimental, or survey research designs⁶. This involves careful selection of independent variables, which might be either manipulated or left untreated, and of dependent variables that act as surrogates for specific behaviours. Once information about all variables is collected, data is quantitatively analysed using a wide range of variance-based methods (e.g., analysis of variance, regression, structural equation modelling).

Behavioural studies that use a variance research methodology can produce a good picture of the generative mechanisms underpinning behavioural processes if they test hypotheses about those mechanisms. For example, variance studies in BOR have examined the impact of individual differences in cognitive motivation and cognitive style on the

conduct of OR-supported activity (Fasolo & Bana e Costa, 2014; Franco et al., 2016b; Lu et al., 2001). There is also a long tradition of testing the behavioural effects of reconfiguring different aspects of OR-supported activity such as varying model or information displays (Bell & O'Keefe, 1995; Gettinger et al., 2013), and preference elicitation procedures (Cavallo et al., 2019; Hämäläinen & Lahtinen, 2016; Pöyhönen et al., 2001; von Nitzsch & Weber, 1993).

A *process methodology* is used to examine OR-supported activity as a series of events that bring about or lead to some behaviour-related outcome. Specifically, it considers as the unit of analysis an evolving individual or group whose behaviour is led by, or leading, the occurrence of events (Poole, 2004). Process explanations take the form of theoretical narratives that account for how event dynamics lead to a final outcome (Poole, 2007). These narratives are often derived from observation, but it is also possible to use an established narrative (e.g., a theory) to guide observation that further specifies the narrative.

Diverse and eclectic research designs are used to implement a process research methodology. Central to these designs is the task of identifying or reconstructing the process through the analysis of events taking place over time. For example, there is an important stream of BOR studies that examines the process of building models by experts and novices (Tako, 2015; Tako & Robinson, 2010; Waisel et al., 2008; Willemain, 1995; Willemain & Powell, 2007). There is also an increasing interest to use process methodologies to take a closer look at actual behaviour in OR-supported settings both, before, during and after OR-related activity is undertaken (Franco & Greiffenhagen, 2018; Käki et al., 2019; Velez-Castiblanco et al., 2016; White et al., 2016).

The variance and process approaches may seem opposite to each other, but instead they should be seen as complementary (Franco et al., 2021; Van de Ven & Poole, 2005). BOR studies using a variance research methodology can explore and test the mechanisms that drive process explanations of behaviour, while BOR studies adopting a process research methodology can explore and test the narratives that ground variance explanations of behaviour. One way of combining a variance and process approach within a single BOR study is by adopting *modelling* as a research methodology. A modelling approach would create models that capture the mechanisms that generate a process of interest such as, for example, trust on an OR-derived solution, and the model can be run to generate the characteristics of that process. Model parameters and structure can then be varied systematically to enable variance-based comparisons
<figure><img src="image_21.png" /></figure>

of trust levels. Furthermore, the trajectory of trust levels over time can be used to gain insights into the nature of the trust development process. As already mentioned, there is a long behavioural modelling tradition within OR but, as far as we know, its potential as a research methodology tool to specifically examine behaviour in OR-supported settings is yet to be realised.

In sum, the variance, process and modelling methodologies offer rich possibilities for the study of human behaviour in OR-supported settings. Which is best for a particular study will depend on the types of question being addressed by BOR researchers, their assumptions about human behaviour, and the data they have access to. Ultimately, a thorough understanding of behaviour in OR-supported settings is likely to require all three research methodologies.

For a detailed review of BOR studies the reader is referred to Franco et al. (2021). A review of behavioural studies in the context of OR in health has been written by Kunc et al. (2020). There are also two collections edited by Kunc et al. (2016) and White et al. (2020). The *European Journal of Operational Research* published a feature cluster on BOR edited by Franco and Hämäläinen (2016a). Finally, BOR-related news and events can be found on the sites of the European Working Group on Behavioural OR<sup>7</sup>, and the UK BOR Special Interest Group<sup>8</sup>.

## 2.3. Business analytics<sup>9</sup>

Business Analytics has its origins in practice, rather than theory, as illustrated by some of the earliest publications on the subject (e.g., Kohavi et al., 2002). Senior executives began to realise the importance of analytics in the first decade of the new millennium because of the ready availability of large amounts of data, the maturity of business performance management, the emergence of self-service analytics and business intelligence, and the declining cost of computing power, data storage and bandwidth (Acito & Khatri, 2014).

Davenport and Harris (2007) gave examples of companies becoming ‘analytical competitors’ by using analytics to support distinctive organisational capabilities. To achieve this level of maturity, it was argued that analytics needs to become a strategic competency. In the 1990s, Fildes and Ranyard (1997) reported on the closure or dispersal of Operational Research groups. Davenport et al. (2010) reflected a reversal of that trend, by focusing on how analytical talent can be organised as an internal resource. They suggested that there are four categories of people to be considered when finding, developing and managing analysts: champions, professionals, semi-professionals and amateurs. In

2012/13, the Institute for Management Science and Operations Research (INFORMS) introduced the Certified Analytics Professional program and examination. This covers the broad spectrum of skills required of analytics professionals, including business problem framing, analytics problem framing, data (handling), methodology selection, model building, deployment and lifecycle management (INFORMS, 2022).

The development of talent is just one of the prerequisites for Business Analytics to create value. Vidgen et al. (2017) recommended ‘coevolutionary change’, aligning their analytics strategy with their strategies for Information and Communications Technology, human resources and the whole business. This helps to ensure that the necessary data assets are available, the right culture is developed to build data and analytics skills, and that there is alignment with the business strategy for value creation. Hindle and Vidgen (2018) proposed a Business Analytics Methodology based on four activities, namely problem situation structuring, business model mapping, business analytics leverage and analytics implementation. They advocated a soft OR approach, Soft Systems Methodology (Checkland & Poulter, 2006), to support structuring and mapping activities.

Many definitions of Business Analytics have been proposed; for a review of early definitions, see Holsapple et al. (2014). According to Davenport (2013), “By analytics we mean the extensive use of data, statistical and quantitative analysis, explanatory and predictive models, and fact-based management to drive decisions and actions” (p. 7). Mortenson et al. (2015) suggested that analytics is at the intersection of quantitative methods, technologies and decision making. Rose (2016) considered analytics as the union of Data Science (which is data centric) and Operational Research (which is problem centric). Power et al. (2018) proposed the following definition: “Business Analytics is a systematic thinking process that applies qualitative, quantitative and statistical computational tools and methods to analyse data, gain insights, inform and support decision-making”. Delen and Ram (2018) pointed out that, although analytics includes analysis, it also involves synthesis and subsequent implementation. These broad perspectives, emphasising synthesis as well as analysis, and qualitative as well as quantitative approaches, are consistent with earlier writings on the use of a broad range of methods in Management Science (e.g., Mingers & Brocklesby, 1997; Pidd, 2009).

Business Analytics can be viewed from different orientations. From a methodological viewpoint, the subject covers descriptive, predictive and prescriptive methods (Lustig et al., 2010). These three
<figure><img src="image_22.png" /></figure>

categories are sometimes extended to four, with a distinction being drawn between ‘descriptive’ and ‘diagnostic’ analytics, following the Gartner analytics ascendancy model (Maoz, 2013). Lepenioti et al. (2020) argue that it is preferable to maintain the threefold categorisation to ensure consistency, with each category addressing both ‘What?’ and ‘Why?’ questions. (Descriptive: ‘What happened?’, ‘Why did it happen?’; Predictive: ‘What will happen?’, ‘Why will it happen?’. Prescriptive: ‘What should I do to make it happen?’, ‘Why should I make it happen?’). For detailed literature reviews on descriptive, predictive and prescriptive analytics, the reader is directed to Duan and Xiong (2015), Lu et al. (2017), and Lepenioti et al. (2020), respectively.

From a technological viewpoint, Business Analytics is facilitated by the integration of transactional data with big data streaming from social media platforms and the Internet of Things into a unified analytics system (Shi & Wang, 2018). These authors suggest that this integration can be achieved in two stages, starting with integration of traditional Enterprise Resource Planning (ERP) and big data, and proceeding to integration of big-data ERP with Business Analytics. Ruivo et al. (2020) reported that analytics ranked second in extended ERP capabilities (behind collaboration) according to the views of 20 experts engaged in a Delphi study. Romero and Abad (2022) suggested that cloud-based big data analytics software will not provide competitive advantage to firms that have not installed a large ERP system, although it will ensure that they do not lag further behind their sector-leading competitors.

From an ethical viewpoint, Business Analytics faces a number of challenges. Davenport et al. (2010) recognised that issues of data privacy can be difficult to address, especially if an organisation operates in a wide range of territories or industries. Ram Mohan Rao et al. (2018) summarised major privacy threats in data analytics, namely surveillance, disclosure, discrimination, and personal embarrassment and abuse, and reviewed privacy preservation methods, including randomisation and cryptographic techniques. A further ethical issue is that AI algorithms are likely to replicate and reinforce existing social biases (O'Neil, 2016). Such algorithmic bias is said to occur when the outputs of an algorithm benefit or disadvantage certain individuals or groups more than others without a justified reason. Kordzadeh and Ghasemaghaei (2022) reviewed the literature on algorithmic bias and showed that most studies had examined the issue from a conceptual standpoint, with only a limited number of empirical studies. Similarly, Vidgen et al. (2020) reviewed papers on ethics in Business Analytics and found that most were at the level of guiding principles and frameworks, with little of

direct applicability for the practitioner. Their case study demonstrated how ethical principles (utility, rights, justice, common good and virtue) can be embedded in analytics development. For further discussions on ethics and OR, the reader is referred to Ormerod and Ulrich (2013), Le Menestrel and Van Wassenhove (2004), and Mingers (2011a) but also §3.8.

Analytics maturity models have been developed to describe, explain and evaluate the development of analytics in an organisation. Król and Zdonek (2020) reviewed 11 maturity models and assessed them in terms of the number of assessment dimensions, scoring mechanism, number of maturity levels, and the public availability of the methodology. They found that the most common assessment dimensions were technical infrastructure, analytics culture and human resources, including staff's analytics competencies. Lismont et al. (2017) undertook a survey of companies, based on the DELTA maturity model (Davenport et al., 2010) of data, enterprise, leadership, targets and analysts. They identified four analytics maturity levels from their survey. The most advanced companies tended to use a wider variety of analytics techniques and applications, to organise analytics more holistically, and to have a more mature data governance policy.

A crucial empirical question is whether Business Analytics adds value to an organisation. An early study on the effect of Business Analytics on supply chain performance was conducted by Trkman et al. (2010). They examined over 300 companies, showing a statistically significant relationship between self-assessed analytical capabilities and performance. Oesterreich et al. (2022) conducted a meta-analysis of 125 firm-level studies, spanning ten years of research in 26 countries. They found evidence of Business Analytics having a positive impact on operational, financial and market performance. They also found that human resources, management capabilities and organisational culture were major determinants of value creation, whereas technological factors were less important.

## 2.4. Combinatorial optimisation¹⁰

A Combinatorial Optimisation (CO) problem consists of searching for the optimal element in a finite collection of elements. More formally, given a set of elements and a family of its subsets, each defining a feasible solution and having an associated value, a CO problem is to find a subset having the minimum (or, alternatively, the maximum) value. The subsets may be proper, like, e.g., in the knapsack problem, or represented by permutations, like, e.g., in the assignment problem (see below). Typically, the feasible solutions are not explicitly listed, but are
<figure><img src="image_23.png" /></figure>

described in a concise manner (like a set of equalities and inequalities, or a graph structure) and their number is huge, so scanning all feasible solutions to select the optimal one is intractable. A CO problem can usually be modelled as an *Integer Program* (IP, see also §2.15) with linear or nonlinear objective function and constraints, in which the variables can take a finite number of integer values.

Consider for example the problem of assigning *n* tasks to *n* agents, by knowing the time that each agent needs to complete each task, with the objective of finding a solution that minimises the overall time needed to complete all tasks (*Assignment Problem*, AP). The solution could obviously be found by enumerating all permutations of the integers 1, 2, ..., *n* and selecting the best one. However, this number is so huge that such approach is ruled out even for small-size problem instances: for *n* = 30, we have *n!* ≈ 2.6 · 10³⁰, and the fastest supercomputer on hearth would need millions of years to scan all solutions. The challenge is thus to find more efficient methods. For example, one of the most famous CO algorithms (the *Hungarian algorithm*) can solve assignment problem instances with millions of variables in few seconds on a standard PC.

The algorithm mentioned above can be implemented so as to solve any AP instance in a time of order $n^3$, i.e., in a time bounded by a polynomial function of the input size. Unfortunately, only for relatively few CO problems we know algorithms with such property, while for most of them (NP-hard problems) the best known algorithms can take, in the worst case, a time that grows exponentially in the size of the instance. In addition, Complexity theory (see also §2.5) suggests that the existence of polynomial-time algorithms for such problems is unlikely. On the other hand, CO problems arise in many industrial sectors (manufacturing, crew scheduling, telecommunication, distribution, to mention a few) and hence there is the prominent and practical need to obtain good quality solutions, especially to large-size instances, in reasonable times.

### 2.4.1. Origins

Many problems arising on graphs and networks (see §2.12) belong to CO (the AP discussed above can be described as that of finding a minimum weight perfect matching in a bipartite graph), and hence the origins of CO date back to the eighteen century. In the following, we narrow our focus to modern CO (see Schrijver, 2005). Its roots can be found in the first decades of the past century, when Central European mathematicians developed seminal studies on *matching problems* (König, 1916), *paths* Menger (1927), and *Shortest Spanning Trees* (SST) (Jarník, 1930; Borůvka, 1926, results independently

rediscovered by Prim, 1957 and Kruskal, 1957). The Fifties produced major results on the AP (Kuhn, 1955; 1956, on the basis of the results by König, 1916 and Egerváry, 1931, also see Martello, 2010), the *Travelling Salesman Problem* (Dantzig et al., 1954), and *Network Flows* (Ford & Fulkerson, 1962), as well as fundamental studies on basic methodologies: *dynamic programming* (DP; Bellman, 1957, see §2.9), *cutting planes* (Gomory, 1958, see §2.15), and *branch-and-bound* (Land & Doig, 1960).

### 2.4.2. Problems and complexity

The most important CO problems, for which we know there are polynomial algorithms, are the basic graph-theory problems mentioned in the previous section. Other important problems, which are relevant both from the theoretical point of view and from that of real-world applications, are instead NP-hard. The main NP-hard CO problems arise in the following areas.

**Scheduling.** Given a set of tasks which must be processed on a set of processors, a scheduling problem asks to find a processing schedule that satisfies prescribed conditions and minimises (or maximises) an objective function, frequently related to the time needed to complete all tasks. This huge area, that includes literally hundreds of problems and variants (mostly NP-hard), is also discussed in §3.27.

**Travelling Salesman Problem (TSP).** Given a weighted (directed or undirected) graph, the problem is to find a circuit that visits each vertex exactly once (Hamiltonian tour) and has minimum total weight. This is one of the most intensively studied problems of CO, and is treated in detail in §2.12.

**Vehicle Routing Problems (VRP).** A VRP is a generalisation of the TSP which consists of finding a set of routes for a fleet of vehicles, based at one or more depots, to deliver goods to a given set of customers by satisfying a set of conditions and minimising the overall transportation cost.

**Facility Location.** These problems require to find the best placement of facilities on the basis of geographical demands, installation costs, and transportation costs, so as to satisfy a set of conditions and to minimise the total cost (see §3.13 for a detailed treatment).

**Steiner Trees.** Given a weighted graph and a subset S of vertices, it is requested to find an SST connecting all vertices in S (possibly containing additional vertices). These problems, which generalise both the shortest path problem and the SST, are treated in detail in §2.12.

**Set Covering.** Given a set of elements and a collection of its subsets, each having a cost, we want to find the least cost sub-collection whose union includes (covers) all the elements.
<figure><img src="image_24.png" /></figure>

**Maximum Clique (MC).** A clique is a complete subgraph of a graph (i.e., it is defined by a subset of vertices all adjacent to each other). Given a graph, the MC problem is to find a clique of maximum cardinality (or, if the graph is weighted, a clique of maximum weight). We refer the reader to §2.12 for a detailed analysis.

**Cutting and Packing (C&P).** Given a set of "small" items, and a set of "large" containers, a problem in this area asks for an optimal arrangement of the items into the containers. Items and containers can be in one dimension (Knapsack Problems (KP), Bin Packing problems) or in more usually two or three - dimensions (C&P). See §3.3 for more details.

**Quadratic Variants of CO problems.** A currently hot research area concerns CO problems whose "normal" linear objective function is replaced by a quadratic one. This greatly increases difficulty: in most cases problems which, in their linear formulation, can be solved in polynomial time (e.g., the AP) or in pseudo-polynomial time (e.g., the KP) become strongly NP-hard.

### 2.4.3. Exact methods for NP-hard problems

For heuristic and approximation algorithms, we refer the reader to §2.13 and §2.5. With the exception of DP methods (§2.9), most exact algorithms for NP-hard CO problems, as well as most commonly used ILP solvers, are based on implicit enumeration. In the worst case, they can require the evaluation of all feasible solutions, and hence computing times growing exponentially with the problem size. The most common methods can be classified as

* Branch-and-Bound (B&B);

* Branch-and-Cut (B&C);

* Branch-and-Price (B&P).

We will describe B&B, the other methods (and their combinations, as B&C-and-Price) being its extensions described in §2.15.

We consider a *maximisation* CO problem having an IP model with inequality constraints of '$\le$' type. For a problem *P*, having feasible solution set *F*(*P*), *z*(*P*) denotes its optimal solution value, and *ub*(*P*) an upper bound on *z*(*P*). The main ingredients of B&B are branching scheme and upper bound computation.

**Branching scheme.** The solution is obtained as follows:

i. subdivide *P* into *m* subproblems, each having the same objective function as *P* and a feasible solution set contained in *F*(*P*), such that the

union of their feasible solution sets is *F*(*P*). The optimal solution of *P* is thus given by the optimal solution of the subproblem having the maximum objective function value;

ii. iteratively, if a subproblem cannot be immediately solved, subdivide it into additional subproblems.

The resulting method can be represented through a *branch-decision tree*, where the root node corresponds to *P* and each node corresponds to a subproblem.

A node of the tree can be eliminated if the feasible solution set of the corresponding subproblem is empty, or its upper bound is not greater than the value of the best feasible solution to *P* found so far.

**Upper bound computation.** A valid upper bound *ub*(*P*) can be computed as the optimal solution value of a *Relaxation* of the IP model of *P*, defined by:

i. a feasible solutions set containing *F*(*P*);

ii. an objective function whose value is not smaller than that of *P* for any solution in *F*(*P*).

A relaxation is "good" if the resulting upper bound *ub*(*P*) is "close" to *z*(*P*) (i.e., if the gap between the two values, *z*(*P*) − *ub*(*P*), is "small"), and the relaxed problem is "easy" to solve, i.e., its optimal solution can be obtained with a computational effort much smaller than that required to solve *P*.

### 2.4.4. Relaxations

The most common relaxation methods are:

* **Constraint elimination:** a subset of constraints is removed from the IP model of *P*, so that the resulting problem is easy to solve. The most widely used case is the linear relaxation;

* **Linear relaxation:** when the model is an *Integer Linear Problem* (ILP), removing the constraints that impose integrality of the variables leads to a *Linear Program* (LP), which is polynomially solvable, commonly used in ILP solvers (see §2.15);

* **Surrogate relaxation:** a subset $\Sigma$ of inequality constraints is replaced by a single **surrogate** inequality, so that the corresponding relaxed problem is easy to solve. The **surrogate** inequality is obtained by multiplying both sides of each inequality of $\Sigma$ by a non-negative constant, and summing, respectively, the left-hand and right-hand sides of the resulting inequalities;

* **Lagrangian relaxation:** a subset $\Lambda$ of inequality constraints is removed from the model and "embedded", in a Lagrangian fashion, into the objective function. For each inequality of $\Lambda$, the difference between left-hand and right-hand sides
(slack) multiplied by a non-negative constant is added to the objective function.

The relaxations can be strengthened by adding one or more valid inequalities (cuts) to the IP model of P, such that they are redundant for the IP model, but can become active when the IP model is relaxed (see §2.15).

### 2.4.5. Further readings

We refer the reader to the following selection of references for more details on the topics covered in this section. Well known, pre-1990 books are those by Garfinkel and Nemhauser (1972, IP), Christofides (1975, algorithmic graph theory), Garey and Johnson (1979, complexity), Burkard and Derigs (1980, AP), Lawler et al. (1985, TSP), and the CO specific volumes by Lawler (1976), Christofides et al. (1979), Papadimitriou and Steiglitz (1982), Martello et al. (1987), and Nemhauser and Wolsey (1988). We list more recent contributions in the order in which the topics were introduced:

* CO: Cook et al. (1998), Schrijver (2003);

* AP: Burkard et al. (2012) for linear and quadratic AP, Cela (2013) for quadratic AP;

* Network Flows: Ahuja et al. (1993);

* Scheduling: Błażewicz et al. (2001, 2007), Pinedo (2012);

* TSP: Gutin and Punnen (2006), Applegate, et al. (2007), Cook (2011);

* VRP: Golden et al. (2008), Toth and Vigo (2014);

* Facility Location: Mirchandani and Francis (1990), Laporte et al. (2015);

* Steiner trees: Hwang et al. (1992), Prömel and Steger (2012). Also see the recent survey by Ljubić (2021);

* Cutting and packing: Martello and Toth (1990), Kellerer et al. (2004). Also see the recent survey by Cacchiani et al. (2022a, 2022b).

## 2.5. Computational complexity¹¹

Operational Research develops models and solution methods for problems arising from practical decision making scenarios. Often, these solution methods are algorithms. The difficulty of a problem can be assessed empirically by evaluating the running times of corresponding algorithms, which requires careful implementations and meaningful test data. Moreover, this can be time-consuming and yields insights that depend on the skills of the programmer and are limited to the available test instances. Computational complexity represents an alternative approach that allows for a more general assessment

of a problem's difficulty that is independent of specific problem instances or solution algorithms.

### 2.5.1. Problem encoding and running times of algorithms

In complexity theory, the running time of an algorithm is expressed in terms of the size of the input, i.e., the amount of data necessary to encode an instance of the problem. Since computers store data in the form of binary digits (bits), the standard binary encoding represents all data of a problem instance in the form of binary numbers. The number of required bits (the encoding length) of an integer is roughly given by the binary logarithm of its absolute value. As an example, consider the binary encoding of instances of the well-known 0-1 knapsack problem (KP). An instance of KP consists of $n$ items – each with a non-negative, integer weight and profit – and a positive, integer knapsack capacity $c$. We can assume that all $n$ item weights are bounded by the capacity $c$ and denote the value of the largest item profit by $p_{max}$. Then, the encoding length of a KP instance is bounded by $(n+1) \cdot \log_2(c) + n \cdot \log_2(p_{max}) \le (2n+1) \cdot \log_2(\max\{c, p_{max}\})$.

Rational numbers can be straightforwardly represented by their (integer) numerator and denominator, but their presence in the input might already influence a problem's computational complexity (Wojtczak, 2018). Irrational numbers cannot be encoded in binary without rounding them appropriately, which means that a different kind of complexity theory is required when general real numbers are part of the input (see Blum et al., 1998, for details). Hence, the following exposition is restricted to the case of integer inputs, where the encoding length of an instance can be bounded by the number of integers needed to represent it multiplied with the binary logarithm of the largest among their absolute values (see the bound for KP instances provided above as an example).

To allow universal running time analyses of algorithms that are independent of specific computer architectures, asymptotic running time bounds described using the so-called $O$-notation (Cormen et al., 2009) are used. Informally, every polynomial in $n$ with largest exponent $k$ is in $O(n^k)$. All terms with exponents smaller than $k$ and the constant coefficient of $n^k$ are ignored. One is then often interested in polynomial-time algorithms whose running time is in $O(|I|^k)$ for some constant $k$, where $|I|$ denotes the encoding length of instance $I$. A less preferred outcome would be a pseudopolynomial-time algorithm, where the running time is only required to be polynomial in the number of integers in the input and the largest among their absolute values (or, equivalently, in the exponentially larger
<figure><img src="image_25.png" /></figure>

encoding length of the input when using *unary encoding*, where the encoding length of an integer is roughly its absolute value).

### 2.5.2. The complexity classes $\mathcal{P}$ and $\mathcal{NP}$

Most application scenarios encountered in Operational Research finally lead to an optimisation problem (often a combinatorial problem – see §2.4), where a feasible solution is sought that minimises or maximises a given objective function. Every optimisation problem immediately yields an associated decision problem, asking a yes-no question. For example, a minimisation problem consisting of a set $\mathcal{X}$ of feasible solutions and an objective function $f$ can be written as $\min\{f(x) : x \in \mathcal{X}\}$. For a given target value $v$, the associated decision problem then asks: Does there exist a feasible solution $x \in \mathcal{X}$ such that $f(x) \le v$? Solving an optimisation problem to optimality trivially answers the associated decision problem for any given $v$. On the other hand, every algorithm for the decision problem can be used to solve the underlying optimisation problem. Given upper and lower bounds, the optimal solution value can be identified in polynomial time by performing binary search between these bounds using the decision problem to answer the query in every iteration of the binary search (assuming that the range of objective function values and the encoding lengths of the bounds are polynomially bounded).

Motivated by the above, the computational complexity of an optimisation problem follows from the complexity of its associated decision problem. Here, the most relevant complexity classes in Operational Research are probably $\mathcal{P}$ and $\mathcal{NP}$, which are often used to draw the line between “easy” and “hard” problems in this context. Formally, the class $\mathcal{P}$ (“polynomial”) consists of all decision problems for which a polynomial-time solution algorithm exists on a deterministic Turing machine (or, equivalently, in any other “reasonable” deterministic model of computation), while the class $\mathcal{NP}$ (“nondeterministic polynomial”) consists of all decision problems for which the same holds on a *nondeterministic* Turing machine. Equivalently, $\mathcal{NP}$ is the class of all decision problems such that, for any yes instance $I$, there exists a certificate with encoding length polynomial in $|I|$ and a deterministic algorithm that, given the certificate, can verify in polynomial time that the instance is indeed a yes instance. Since the most natural certificate is often a (sufficiently good) solution of the problem, $\mathcal{NP}$ can informally be defined as the class of decision problems for which solutions can be verified in polynomial time. For example, when considering the travelling salesman problem (TSP) on a given edge-weighted graph, the associated decision problem asks whether or not there exists a tour

(Hamiltonian cycle) of at most a given length $v$. While no polynomial-time algorithm for this decision problem is known to date, the problem can easily be seen to be in $\mathcal{NP}$ since the natural certificate to provide for a yes instance is simply a tour with length at most $v$, whose feasibility and length can be easily verified in polynomial time.

Observe that these definitions directly imply that $\mathcal{P} \subseteq \mathcal{NP}$. Most researchers believe that $\mathcal{P} \subseteq \mathcal{NP}$ or, equivalently, that there are problems in $\mathcal{NP}$ that do not admit polynomial-time solution algorithms. However, formally proving that $\mathcal{P} \ne \mathcal{NP}$ (or that $\mathcal{P} = \mathcal{NP}$) is still one of the most famous open problems in theoretical computer science to date.

This so-called *P versus NP problem* can be equivalently expressed using the well-known notion of *NP-completeness* (see, e.g., Garey & Johnson, 1979). Intuitively, *NP-complete* problems are the hardest problems in *NP* in the sense that, if one of these problems admits a polynomial-time solution algorithm, then so does every problem in *NP* (and, thus, we would obtain $\mathcal{P} = \mathcal{NP}$). A decision problem (not necessarily in *NP*) with this property is also called *NP-hard*. This means that a problem is *NP-complete* if and only if it is both *NP-hard* and contained in *NP*. The first problem shown to be *NP-complete* in Cook's famous theorem (Cook, 1971) is the (Boolean) satisfiability problem (SAT). Shortly after, Karp (1972a) gave a list of 21 fundamental problems that are *NP-complete*. While Cook's proof that SAT is *NP-complete* required considerable effort, proving that further problems are *NP-complete* became significantly easier with this knowledge and hundreds – if not thousands – of problems were shown to be *NP-complete*.

A decision problem is *NP-complete* if and only if (1) it is contained in $\mathcal{NP}$ and (2) some $\mathcal{NP}$-complete problem (and, therefore, all problems in $\mathcal{NP}$) can be reduced to it via a *polynomial-time reduction*. Such a polynomial-time reduction works as follows: For any instance of the known $\mathcal{NP}$-complete problem (e.g., SAT or TSP), one has to construct an instance of the investigated problem in polynomial time such that the two instances are equivalent, i.e., the constructed instance is a yes instance if and only if the given instance is a yes instance. Note that the requirement that the instance must be constructed in polynomial time (and, therefore, have encoding length polynomial in the encoding length of the original instance) is crucial. A common error in reductions is that the encoding length of the constructed instance depends polynomially on the size of numerical values in the given instance (instead of their encoding length).

The importance of the encoding can be illustrated by the 0-1 knapsack problem (KP), which is *NP-hard* if binary encoding is used, but can be
<figure><img src="image_26.png" /></figure>

solved in polynomial time (via dynamic programming) if unary encoding is used (so $\mathcal{NP}$-hardness of the unary-encoded version would imply that $P = \mathcal{NP}$). Problems like this, i.e., problems whose binary-encoded version is $\mathcal{NP}$-hard, but whose unary-encoded version can be solved in polynomial time, are called weakly $\mathcal{NP}$-hard, while problems (such as SAT) that remain $\mathcal{NP}$-hard also under unary encoding are called strongly $\mathcal{NP}$-hard. The existence of a pseudopolynomial-time algorithm is possible for weakly $\mathcal{NP}$-hard problems, but not for strongly $\mathcal{NP}$-hard problems (unless $P = \mathcal{NP}$).

### 2.5.3. Approximation algorithms

While some realistic-size instances of $\mathcal{NP}$-hard problems might still be solvable in reasonable time, this is not the case for all instances. In general, one can deal with $\mathcal{NP}$-hardness by relaxing the requirement of finding an optimal solution and instead settling for a “good-enough” solution. This leads to heuristics, whose aim is producing good-enough solutions in reasonable time (see §2.13 for details) and approximation algorithms (Vazirani, 2001; Williamson & Shmoys, 2011; Ausiello et al., 1999). Given $\alpha \ge 1$, an $\alpha$-approximation algorithm for an optimisation problem is a polynomial-time algorithm that, for each instance of the problem, produces a solution whose objective value is at most a factor $\alpha$ worse than the optimal objective value. The factor $\alpha$, which can be a constant or a function of the instance size, is then called the approximation ratio or performance guarantee of the approximation algorithm. While it is standard to use $\alpha \ge 1$ for minimisation problems, there is no clear consensus in the literature as to whether $\alpha \ge 1$ or $\alpha \le 1$ should be used for maximisation problems. For example, the simple extended greedy algorithm for the knapsack problem produces a solution with at least half of the optimal objective value on each instance, i.e., it is a 1/2- or a 2-approximation algorithm.

While inapproximability results can be shown for some $\mathcal{NP}$-hard problems (see Hochbaum, 1997, ch. 10), others allow for approximation algorithms with approximation ratios arbitrarily close to one, i.e., they admit a polynomial-time approximation scheme (PTAS). A PTAS is a family of algorithms that contains a $(1 + \varepsilon)$-approximation algorithm for every $\varepsilon > 0$. If the running time is additionally polynomial in $1/\varepsilon$, the PTAS is called a fully polynomial-time approximation scheme (FPTAS). If all objective function values are integers, every FPTAS can be turned into a pseudopolynomial-time exact algorithm, so strongly $\mathcal{NP}$-hard problems do not admit an FPTAS (unless $P = \mathcal{NP}$). Conversely, pseudopolynomial-time algorithms, in particular dynamic programming

algorithms, often serve as a starting point for designing an FPTAS (Woeginger, 2000; Pruhs & Woeginger, 2007).

### 2.5.4. Further complexity classes

Theoretical computer science developed a wide range of complexity classes far beyond the $P$ vs. $\mathcal{NP}$ dichotomy. Considering algorithms requiring polynomial space, i.e., for which the encoding length of the data stored at any time during the algorithm's execution is polynomial in the encoding length of the input (but no bound on the running time is required), gives rise to the class PSPACE. It is widely believed that $\mathcal{NP} \subsetneq PSPACE$, but even whether $P \neq PSPACE$ holds is not known.

In the theoretical analysis of bilevel optimisation problems (see, e.g., Labbé & Violin, 2016) the complexity class $\Sigma_2^p$ plays an important role (see Woeginger, 2021). Here, a yes instance $I$ is characterised by the existence of a certificate of encoding length polynomial in $|I|$, such that a certain polynomial-time-verifiable property holds true for all elements of a given set $\mathcal{Y}$. As an example, consider the 2-quantified (Boolean) satisfiability problem. Here, an instance consists of two sets $X$ and $Y$ of Boolean variables and a Boolean formula over $X \cup Y$. The question then is whether there exists a truth assignment of the variables in $X$ such that the formula evaluates to true for all possible truth assignments of the variables in $Y$. This definition immediately sets the stage for a bilevel problem, where the decision $x$ of the upper level (the leader) should guarantee a certain outcome for every possible decision $y$ at the lower level (the follower). It is widely believed that $\mathcal{NP} \subsetneq \Sigma_2^p$ although $\Sigma_2^p$-hardness does not rule out the existence of a PTAS (Caprara et al., 2014). Under this assumption, $\Sigma_2^p$-hardness does, however, rule out the existence of a compact ILP-formulation, which can be a valuable finding for bilevel optimisation problems.

For some $\mathcal{NP}$-hard problems, one can construct algorithms with running time $O(f(k) \cdot \text{poly}(|I|))$ for an arbitrary computable function $f$, where the parameter $k$ describes a property of the instance $I$. Such problems are called fixed-parameter tractable. For example, the satisfiability problem SAT is fixed-parameter tractable with respect to the parameter $k$ that represents the tree-width of the primal graph of the SAT instance. In this graph, the vertices are the variables and two vertices are joined by an edge if the associated variables occur together in at least one clause, see Szeider (2003). This parametric point of view is captured in the W-hierarchy of complexity classes – see Niedermeier (2006) and the seminal book by Downey and Fellows (1999).
<figure><img src="image_27.png" /></figure>

## 2.6. Control theory¹²

Control theory deals with designing a control signal so that the state or output variables of the system meet certain criteria. It is a broad umbrella term that covers a variety of theories and techniques. Control theory has been widely applied in the studies of economics (Tustin, 1953; Grubbström, 1967), operations management (Simon, 1952; Vassian, 1955, also see (Sarimveis et al., 2008) for a recent review), and finance (Sethi & Thompson, 1970). Here, we do not intend to provide an exhaustive or comprehensive review. Instead, we try to structurally organise the concepts and techniques commonly applied in operations research, which means that technical details will be omitted. We direct interested readers to a number of textbooks in the reference list, and an excellent review by Åström and Kumar (2014) for those interested in the development of control theory.

The major distinction between control theory and other optimisation theories is that the control variable to be designed is normally a time-varying, dynamic function. The control signal can either be dependent on the state variables (which is referred to as feedback control or closed-loop control) or independent (feedforward control or open-loop control). The design of control signals and control policies (defined as the function between the state of the system and the control, also known as “control laws” or “decision rules”) is based on the structure of the system to be controlled (sometimes called the “plant” in the control engineering literature). Thus, the type of the dynamical system often define the type of control problem. In continuous systems, the time variable is defined on the real axis, suitable to describe continuous processes such as fluid processing and finance. In discrete systems, time is defined as integers, suitable in cases such as production and inventory control, where the production quantity is released every day. Linear systems are comprised of linear (or affine) state equations, while nonlinear systems contain nonlinear elements. Nonlinear systems are more difficult to analyse and control, and may lead to complex system behaviours such as bifurcation, chaos and fractals (Strogatz, 2018). But there are linearisation strategies which approximate the nonlinear system locally as linear systems (Slotine et al., 1991). Based on whether random input is present, the dynamical system can be categorised into deterministic and stochastic.

There are two fundamental methods in the analysis of the system and control. The first relies on time-frequency transformations (Laplace transform for continuous systems and z-transform for discrete systems). A transfer function in the frequency domain can be used to represent and analyse the system (Towill, 1970). This method saves

computational effort; however, it can only deal with linear system models and each transfer function only describes the relation between a single input and a single output (SISO). The second method directly tackles the state equations in the time domain and describes the movement of system state in the state space. It is suitable for nonlinear systems and multi-input-multi-output (MIMO) systems. With the advancement of computing technology, the computational burden faced by the time-based method becomes less significant. The literature refers to the frequency-based method as *classic control theory* (Ogata et al., 2010) and the time-based method as *modern control theory*.

The system under the effect of the control policy must be examined with respect to its properties and dynamic performance. *Stability* is a property of the dynamical system, that the system can return to its steady state after receiving a finite external disturbance. Stability is a fundamental precondition that almost all control designs must meet, with few exceptions such as clocks and metronomes, where a periodic or cyclic response is desired. The stability criterion is straightforward to derive for linear systems, where both frequency-based (e.g., Routh-Hurwitz stability criteria and Jury's inners approach, Jury & Paynter, 1975) and time-based (e.g., the eigenvalue approach) methods exist. However, stability analysis for nonlinear systems is more challenging (Bacciotti & Rosier, 2005). Other important properties of the control system include *controllability*, defined as the ability to move the system to preferred state using only the control signal; and *observability*, defined as the ability to infer the system state using the observable output signals (Gilbert, 1963).

In addition to these intrinsic properties, the system can also be evaluated by the system's response to some characteristic input functions. The *step function* (sometimes referred to as the Heaviside function) takes the value of zero before the reference time point, and one thereafter. The *impulse function* (the Dirac δ function) takes the value of infinity at the reference time point and zero otherwise. These two input functions usually represent an abrupt change in the external environment. The *sinusoidal function* can be used to describe the periodic and seasonal externalities. The Bode plot describes the amplitude and phase shift between the sinusoidal input and output. For stochastic environments, the *white noise* is used to mimic random disturbances. It is a random signal that follows an identical and independent Gaussian distribution and has a constant power spectrum. The noise bandwidth of the system determines the ratio between output and input variances when the input is iid. The value of the noise bandwidth can be derived from either the transfer function or the state
<figure><img src="image_28.png" /></figure>

space representation. This concept is used in analysing the amplification phenomena in supply chains (see §3.24).

In practise, the system state and even the system structure may be unknown. Therefore, statistical techniques, known as *state estimation* and *system identification*, have been developed. State estimation uses observable output data to estimate the unobservable system states. A popular technique for this purpose is the Kalman filter (Kalman, 1960), essentially an adaptive estimator that can be applied not only in linear, time-invariant cases (LTI, where the system is linear and does not change over time), but also non-linear and time-variant cases. For example, it has been applied to estimate the demand process from the observed sales data (Aviv, 2003). System identification attempts to “guess” the structure of the system from the input and output.

Along with the development of control theory, various control strategies have been proposed. They are designed to fit the structure of the system, the objective of the control, and most importantly, to offer a paradigm to design the control policy. In what follows, we provide a brief summary of control strategies. *Linear* control strategies can be represented linearly (in the form of transfer function). They offer great analytical tractability and satisfactory performance, especially when the open-loop system is also linear. Two widely adopted policies in this family are proportional-integral-derivative (PID) control and full-state feedback (FSF) control. In PID control, an error signal between the output and the reference input (e.g., a Heaviside function) is computed. The control signal is a linear combination of the error, the integral, and the derivative of the error. These three components can appear separately. The proportional control has been applied in mechanical and managerial mechanisms such as the centrifugal governor and production planning (Chen & Disney, 2007). The full-state feedback control defines the control signal as a linear combination of the full system state vector, where the coefficient vector (the “gain”) shares the same dimensionality as the state vector. By tuning the gain, the poles of the closed-loop system (the eigenvalue of the transition matrix or the roots of the characteristic equation) will change their position in the complex plane, adjusting the system performance. The full-state feedback policy can also be applied in production and inventory control (Gaalman, 2006).

In contrast to the linear strategies, the *nonlinear* control strategies are defined as policies where the control signal cannot be represented by a linear function of the system state (Slotine et al., 1991). These policies are primarily used when the open loop system is also nonlinear. One such policy is sliding mode control, where the control signal is a

switching function of the state, dependent on some switching rules. The system is then maintained near a hypersurface of the state space (sliding), where the dynamic behaviour of the system is desired. It should be ensured that the hypersurface is reachable from any initial state and that the system state can be maintained on the hypersurface by the policy. In practise, bang-bang control is adopted frequently as a special case of sliding mode control, where the control signal can take only two possible values. The rocket engines and domestic thermostats are examples of such (with on and off states).

*Optimal control* aims at finding the control signal or control policy that allows an objective function to reach its extreme point (Sethi & Thompson, 2009; Bertsekas, 2012a). The objective function could be dependent on the state, output and/or control. Many control policies mentioned above, e.g., full-state feedback control and sliding mode control, have been proved to be the optimal control of some control problems. Optimal control in the special sense is based on Pontryagin's Maximum (or equivalently Minimum) Principle and mainly deals with the design of the open-loop control signal. When equipped with the Hamilton-Jacobi-Bellman (HJB) theory, it can be used to design optimal feedback control policies. Optimal control is closely connected with dynamic programming, which will be reviewed in §[dynamic programming]. The optimal control technique has been widely applied in operations management (e.g., Kumar & Swaminathan, 2003).

When random external disturbances are present, the *stochastic control* techniques are necessary (Åström, 2012). In these situations, objective functions are usually statistical functions of the state or the output, such as the absolute mean or variance. The most well-studied stochastic control problem is the Linear Quadratic Gaussian (LQG) problem, where the system is linear, the objective function is of quadratic form, and the noise signal follows a Gaussian distribution. The optimal control policy in this case is a linear one. Many supply chain management problems can be modelled in LQG form (Lee et al., 1997). For more complex problems involving nonlinearity or an unspecified system structure, the *model predictive control* (MPC) approach can be used (Camacho & Bordons, 2013). This approach transforms the infinite-horizon problem into a finite-horizon problem by focusing only on *T* periods in the future, deriving the control signal for these *T* periods, and adopting the most recent control. In the next period, the prediction is updated, and this process is repeated. MPC is not an optimal control method due to the finite-horizon approximation, yet it works very well in practise (Doganis et al., 2008). To deal with parametrical uncertainties in the disturbance, *robust*
<figure><img src="image_29.png" /></figure>

control provides guaranteed performance (Zhou & Doyle, 1998). The well-known $H_\infty$ control (H infinity) is one of such examples. It minimises the largest singular value of the transition matrix function, which in SISO systems equates to the peak value of the frequency response curve. This minimax strategy ensures that any frequency component in the input will not be amplified too much. Finally, if the system parameters vary over time, adaptive control allows the control policy to update according to the estimated parameters (Åström & Wittenmark, 2013). The difference between adaptive and robust control is that the policy is dynamic in the former and static in the latter.

Recent development of control theory can be seen in the controlling of complex, large scale and network system; the use of artificial intelligence in control engineering; and the application of control theory in areas of physics, biology and economics.

## 2.7. Data envelopment analysis¹³

Data Envelopment Analysis (DEA) is a non-parametric frontier analysis methodology mainly used to assess the relative efficiency of a set of homogeneous operating units (termed Decision Making Units, DMUs). DMUs are assumed to consume inputs (i.e., resources) to produce outputs (e.g., goods and services). The production function that indicates the amount of outputs that can be produced from a given input vector is unknown. DEA does not make any assumption about the functional form of that dependency. Instead, DEA uses the observed data to infer the Production Possibility Set (PPS), also called the DEA technology, which contains all the operating points that are deemed feasible. This is achieved on the basis of a few assumptions (like envelopment of the observations, free disposability of inputs and outputs, convexity and returns to scale) and invoking the Minimum Extrapolation Principle. The resulting PPS contains all linear combinations of the observations along with all the operating points that they dominate. This leads to Linear Programming models whose main decision variables are the intensity variables used to compute the target operating point (projection). This target operating point must dominate the DMU being projected and represents maximal improvements (i.e., input reduction and output increase) with respect to the latter. Hence, the computed target belongs to the efficient frontier (which is the non-dominated subset of the PPS) and the efficiency score is a decreasing function of the distance from the DMU to the computed efficient target. There are different ways of measuring this distance, which, ultimately, depends on the potential input and output improvements (i.e., slacks)

computed by DEA. Before diving into the DEA methodology note that, as Cook et al. (2014) point out, although DEA has a strong link with production theory in economics, it is often used to benchmark the performance of manufacturing and service operations. In such benchmarking exercises, the efficient DMUs, as defined by DEA, may not necessarily form a “production frontier”, but rather a “best-practice frontier”. Thus, the purpose of the performance measurement exercise affects the classification of the different variables considered into inputs or outputs.

### 2.7.1. Efficiency assessment and target setting
DEA models

The seminal DEA models by Charnes et al. (1978) and Banker et al. (1984) were oriented (i.e., gave priority to reducing the inputs or to increasing the outputs) and looked for a uniform (i.e., radial) improvement in all the input or output dimensions. The projection can also be estimated using a given direction, giving rise to Directional Distance Function (DDF) DEA models (Wang et al., 2019a). However, most DEA approaches are non-radial and non-oriented (e.g., Fukuyama & Weber, 2009). Actually, because DEA aims at simultaneously improving inputs and outputs, it is inherently a multiobjective optimisation approach. Hence, taking into account the preferences of a decision maker, any Pareto optimal point can be selected as efficient target (Soltani & Lozano, 2020).

Most DEA models compute targets that can be sometimes far away from the observed DMU. This increases the difficult and effort required to achieve the target. Hence, DEA models that compute closest efficient targets have been developed (Aparicio et al., 2007). An alternative is to compute stepwise efficiency improvement approaches that may eventually achieve ambitious efficient targets but after several gradual improvement steps (Lozano & Villa, 2005).

DEA models for handling non-discretionary variables (Banker & Morey, 1986), undesirable outputs (Kuosmanen, 2005), integer variables (Kazemi Matin & Kuosmanen, 2009), ratio variables (Olesen et al., 2022), negative data (Sharp et al., 2007), and fuzzy data (Arana-Jiménez et al., 2022) have also been proposed. Each of the above “complications” requires specific adaptations of the methodology and being capable of taking them into account is a proof of the power and flexibility of DEA.

The DEA models based on the PPS concept are labelled as envelopment formulations. There are also dual multiplier formulations in which the decision variables are not the intensity variables used to compute the target inputs and outputs but the
<figure><img src="image_30.png" /></figure>

corresponding input and output shadow prices. Multiplier formulations let each DMU choose these input and output weights so that its efficiency is maximised. This freedom often leads to DMUs choosing idiosyncratic or unreasonable weights. Imposing Assurance Regions (AR) and other types of weight restrictions has been proposed (Allen et al., 1997) as well as measuring the efficiency of the DMUs as the average of the cross-efficiency scores computed with the input and output weights chosen by the different DMUs (Doyle & Green, 1994; Chen & Wang, 2022). Another alternative that has been proposed is using a Common Set of Weights (CSW) instead of letting each DMU choose its own (Salahi et al., 2021).

In addition to computing efficiency scores, DEA can be used to rank the DMU. The problem here is that in conventional DEA all the DMUs labelled as efficient are tied and cannot be ranked. In addition to the CSW or cross-efficiency approaches mentioned above, there are other DEA-based full ranking methods, like the super-efficiency approach (Tone, 2002). Alternatively, instead of fully ranking the DMUs, ranking intervals and dominance relations can be established (Salo & Punkka, 2011).

## 2.7.2. Dynamic and network DEA models

DEA views DMUs as input-output black boxes. However, it is often the case that DMUs have an internal structure with different stages or processes (sometimes labelled subDMUs). Many different Network DEA (NDEA) models have been developed to address these scenarios (Tone & Tsutsui, 2009). The key features of NDEA models are that each process has its own technology and that, except in the case of parallel processes, there exist intermediate product flows between the processes. Some NDEA models can compute an efficiency score for each process and relate the overall efficiency score to the scores of the individual processes (Kao, 2016). It must be noted that the NDEA configuration most studied and most commonly used in practice involves two stages in series (see, e.g., Cook, et al., 2010; Halkos et al., 2014).

Multi-period and dynamic scenarios can be modelled in a manner similar to NDEA simply by considering each time period as a subDMU. The difference between multi-period approaches (Kao & Liu, 2014) and Dynamic DEA (Tone & Tsutsui, 2010) is that in the latter there are flows between consecutive periods (i.e., carryovers). Dynamic NDEA (DNDEA) models, in which there are carryovers between periods as well as intermediate product flows between the processes, have also been developed (Tone & Tsutsui, 2014).

## 2.7.3. Centralised DEA models

DEA generally projects each DMU separately onto the efficient frontier. There are situations in which the DMUs belong to the same organisation and there is a Central Decision Maker (CDM) that is interested in the overall system performance and therefore in projecting all the DMUs simultaneously. This type of Centralised DEA (CDEA) models are commonly used for resource allocation (Lozano & Villa, 2004) and for centralised production planning (Lozano, 2014). Also, an approach to measure the centralised efficiency of the individual DMUs in CDEA scenarios has been proposed (Davtalab-Olyaie et al., 2023).

DEA models for allocating a fixed input or common revenue (Li et al., 2021) or for fixed-sum-outputs (FSO; Zhu et al., 2017) also share with CDEA the need to project all the DMUs simultaneously to take into account their interrelationships. These models, same as CDEA, can use an envelopment or a multiplier formulation. While the key feature of the former is that all DMUs are projected simultaneously, that of the latter is that, same as in CSW, a single set of input and output weights is considered.

## 2.7.4. DEA and total factor productivity (TFP) growth

DEA can be used to compute the Malmquist Productivity Index (MPI) by projecting the DMU in two consecutive periods onto the efficient frontier of each period and computing the geometric mean of the change in the corresponding radial efficiency scores between the two periods (Färe et al., 1992). The Malmquist-Luenberger Productivity Indicator (MLPI) is analogous but it employs the arithmetic average and an additive decomposition of DDF efficiency scores (Chambers et al., 1996). In both cases, the TFP growth of each DMU can be decomposed into an efficiency change and a technological change component. Other alternative decompositions of the MPI and MLPI have been developed (Epure et al., 2011).

Other approaches compute a global MPI (Pastor & Lovell, 2005; Kao & Liu, 2014). These have the circularity property, missing in the adjacent-periods MPI. Changes in prices can be also incorporated to compute and decompose a global cost MPI (Tohidi et al., 2012). MPI variants that take into account the projections of all the observations or of different groups of observations as well as approaches to compute and decompose the aggregate productivity growth index of a whole industry and input-specific productivity growth indexes have also been proposed (Aparicio et al., 2017; Kapelko et al., 2015).
<figure><img src="image_31.png" /></figure>

### 2.7.5. Metafrontier analysis

In scenarios where the DMUs are heterogeneous and belong to different groups, not necessarily disjoint, the DMUs can be projected onto its group frontier as well as onto the metafrontier that results from enveloping all the group frontiers. Measuring the difference between the corresponding efficiency scores can be used to estimate the distance between both frontiers and hence the corresponding technology gap of each group. Although the group technologies are generally convex, the metatechnology is generally non-convex (Afsharian & Podinovski, 2018).

The metafrontier approach can be used in DNDEA (See et al., 2021) and CDEA (Gan & Lee, 2022) contexts. Also, using metafrontier concepts with each group of observations corresponding to a different time period, meta-MPI and meta-MLPI can be computed and appropriately decomposed (Portela & Thanassoulis, 2010).

### 2.7.6. Other DEA approaches

There are other interesting DEA approaches that have not been covered above, like congestion (Ren et al., 2021), window analysis (Peykani et al., 2021), etc. Moreover, the field, although mature, is still expanding, with promising new developments, like Efficiency Analysis Trees (EAT) (Esteve et al., 2020), Support Vector Frontiers (SVF) (Valero-Carreras et al., 2022), or big data DEA (Dellnitz, 2022). This is not to mention the large and increasing number of DEA applications (see §3.6, §3.7, and §3.19). For further learning on DEA the interested reader is referred to existing textbooks (Cooper et al., 2007), handbooks (Cooper, et al., 2011; Cook & Zhu, 2014; Zhu, 2015) and review papers (Kao, 2014; Contreras, 2020; Peykani, et al., 2020).

## 2.8. Decision analysis¹⁴

The term decision analysis was introduced by Howard (1966) as “a logical procedure for the balancing of the factors that influence a decision”, pointing out that “the procedure incorporates uncertainties, values and preferences in a basic structure that models the decision”. According to Keeney (1982) decision analysis is a “formalisation of common sense for decision problems which are too complex for informal use of common sense” and, in more technical form “a philosophy, articulated by a set of logical axioms, and a methodology and collection of procedures, based upon those axioms, for responsibly analysing the complexities inherent in decision problems”. In a slightly different perspective, Roy (1993) proposed the concept of decision aiding as “the activity of one who, in ways we call scientific, helps to obtain elements of

answers to questions asked by actors involved in a decision-making process, elements helping to clarify this decision in order to provide actors with the most favourable conditions possible for that type of behaviour which will increase coherence between the evolution of the process, on the one hand, and the goals and/or systems of values within which these actors operate on the other”.

For Howard (1966) “the essence of the procedure is the construction of a structural model of the decision in a form suitable for computation and manipulation”. For Keeney (1982) “the foundations of decision analysis are provided by a set of axioms … which provide principles for analysing decision problems”. Moreover, “the philosophical implications of the axioms are that all decisions require subjective judgements and that the likelihoods of various consequences and their desirability should be separately estimated using probabilities and utilities, respectively”. In this perspective, the key components of a decision problem are the set of alternatives to be taken into consideration; the set of consequences describing outcomes of alternatives, possibly in terms of a plurality of attributes or criteria; if the consequences are uncertain, the beliefs about their possible realisations expressed in terms of a probability distribution; the preferences of the decision maker. The objective of the decision analysis is to construct a value function representing the preferences of the decision maker by assigning each alternative an evaluation of its desirability. In case of uncertainty of the consequences, the value function is expressed in terms of expected value with respect to the probability of the consequences. The basic methodology to induce the value function is based on the pioneering work of von Neumann and Morgenstern (1944) that showed that a small set of axioms imply that the “utility” of an outcome $x$ is defined as the probability of getting the most-preferred outcome and otherwise the least-preferred outcome that would be indifferent to receiving outcome $x$ with certainty. For Roy (1993), the decision aiding procedure should be developed in a constructive approach in which “concepts, models, procedures and results are here seen as suitable tools for developing convictions and allowing them to evolve, as well as for communicating with reference to the bases of these convictions”. In this perspective the “object is not to know or to approximate the best possible decision but to develop a corpus of conditions and means on which we can base our decisions in favour of what we believe to be most suitable”.

Decision Analysis is mainly based on concepts and tools related to the subjective probability of Ramsey (1931) and de Finetti (1937), the theory of expected utility of von Neumann and Morgenstern
<figure><img src="image_32.png" /></figure>

(1944) and subjective expected utility of Savage (1954), the Multiple Attribute Utility Theory (MAUT) of Keeney and Raiffa (1976) and the psychology of judgement and decision-making of Tversky and Kahneman (1974). The general idea is to try to evaluate each alternative by assigning a value based on the utilities of the outcomes obtained in each state of the nature multiplied by their probabilities. Delayed consequences may be discounted according to the time at which they are obtained. Each outcome may be evaluated by considering value trade-offs among multiple attributes. Decision analysis techniques include Utility Function Elicitation techniques, Probability Elicitation protocols, Net Present Value, Decision Trees, Influence Diagrams, and Monte Carlo simulation-based decision analysis (Clemen, 1996); Value-Focused Thinking (Keeney, 1996a); Portfolio Decision Analysis (Salo et al., 2011), Bayesian Networks (Pearl, 1988), and multi-stage decision optimization techniques such as dynamic programming and reinforcement learning.

Considering the distinction between normative, descriptive and prescriptive approaches (Bell et al., 1988), the general perspective of the decision analysis is prescriptive rather than normative or descriptive (Edwards et al., 2007). Descriptive analysis concerns the representation and prediction of observed decisions and normative analysis concerns the decisions that ideally coherent and rational individuals should take. Instead, prescriptive analysis tries to propose methods and techniques that will help real people make better decisions with lower regret and greater coherence of values and behaviors. In this context, decision analysis takes a prescriptive approach that, focusing on the few basic axioms underlying subjective expected utility, adopts “pragmatically” the aspiration to the rationality of the normative approach, trying to correct all the heuristics and biases discovered and investigated by descriptive analysis (Tversky & Kahneman, 1974). The decision aiding approach (Roy, 1993) takes a different perspective that, criticising the idea that there is an objectively optimal decision to be discovered or at least approximated, aims to provide a recommendation consisting in a set of convictions constructed in the course of a decision process based on multiple interactions between the analyst and the decision maker. The decision aiding approach leads directly to a multi-criteria perspective (Belton & Stewart, 2002; Greco et al., 2016) taking explicitly into consideration the multiple attributes or criteria (e.g., related to finance, resources, time, and environmental impacts) to be considered in the decision problem at hand. This avoids the risk of a fictitious, not reasoned and arbitrary conversion of evaluations on different

criteria to a common unit, facilitating the discussion on the respective role of each criterion (Roy, 2005, 1996). To compare alternatives in a multicriteria decision procedure four main approaches can be adopted:

* aggregating criteria assigning a single value to each alternative: this is the case of above mentioned MAUT, as well as of some of the most well known multicriteria methods such as SMART (Edwards & Barron, 1994), and UTA (Jacquet-Lagreze & Siskos, 1982); a specific mention deserves in this context the AHP approach (Saaty, 1977), that is probably the most adopted (although controversial; see, e.g., Dyer, 1990) multicriteria method. It is based on the comparison of “importance” of criteria and of evaluation of alternatives with respect to considered criteria by means of a nine point qualitative scale; another specific class in this family are the distance-based methods which, following the main principle of TOPSIS (Hwang & Yoon, 1981), the first and most famous of these methods, evaluate each alternative on the basis of their distance from the positive ideal solution and the negative ideal solution (the fictitious alternatives that have the best and the worst evaluation on each criterion, respectively); two other well-known methods in this class are VIKOR (Opricovic & Tzeng, 2004) and TODIM (Gomes & Lima, 1991).

* aggregating criteria by means of one or more synthesising preference relations: the most well known methods based on this approach are the ELECTRE methods (Figueira et al., 2013, 2016), that build a crisp or valued preference relation called outranking for which an alternative *a* is at least as good as another alternative *b* if a *a* is not worse than *b* for a majority of important criteria (concordance) and there is no criterion for which the advantage of *b* over *a* is so large that it prevents the possibility to declare *a* at least as good as *b* (non-discordance);

* aggregating criteria through “if ..., then ...” decision rules (Greco et al., 2001): the alternatives obtain an overall evaluation by matching decision rules with a syntax “if the alternative is at least at level $l_{j1}$ on criterion $g_{j1}$ and ... at least at level $l_{jr}$ on criterion $g_{jr}$, then the alternative is globally at least at level $l_{tot}$”, such as “if the student has an evaluation at least good on mathematics and at least medium on literature, then the student is globally at least medium”; these rules are induced from a set of examples of decisions supplied by the decision maker. The advantage of this approach is its explainability due to the fact that the decision rules are expressed in natural language;
<figure><img src="image_33.png" /></figure>

* aggregating criteria through an interactive multi-objective optimisation (Branke et al., 2008): with this approach one can handle decision problems in which a set of objectives have to be optimised under given constraints (see Sawaragi et al., 1985; Steuer, 1985; Miettinen, 1999; Ehrgott, 2005). In this context, the concept of Pareto efficient solution is fundamental: it is a solution for which one cannot improve one objective without deteriorating some others. Several algorithms have been proposed for Pareto set generation and among them let us remember the weighted sum method, the lexicographic method, the achievement secularising function, the epsilon constraint method (for surveys see Marler & Arora, 2004, or Chapters 18 and 19 in Greco et al., 2016). Dealing with a multiobjective optimisation problem, it is important to discover the set of Pareto efficient solutions most preferred by the decision maker. Recently, beyond many exact methods, some heuristic methods have been proposed for these problems, such as some hybridisation between evolutionary multiobjective optimisation algorithms aiming to approximate the whole set of Pareto efficient solutions (Deb, 2001) and some multicriteria preference elicitation methods to guide the optimisation algorithm toward the most interesting set of Pareto efficient solutions (see, e.g., Phelps & Köksalan, 2003; Branke et al., 2016).

## 2.9. Dynamic programming¹⁵

Dynamic programming (DP) was the brainchild of Richard Bellman (Bellman, 1953), who wrote “DP is a mathematical theory devoted to the study of multistage processes”. Indeed, in the seven decades since his seminal work, the uses of DP have grown substantially thanks to its algorithmic nature in solving sequential decision-making problems, where the preceding actions and their realisation (in terms of consequences) will impact on the course of futures. Examples of such problems include multiperiod inventory management, or asset allocation (portfolio management) over a given time horizon. The central idea of DP is to break down the original multistage problem into a number of tail sub-problems by stages. For each stage, the tail sub-problem is a truncated version of the original problem starting from this stage. These tail sub-problems are then recursively solved one by one from the last stage backwards to the first one, at which point the original problem is solved. The solution of such a procedure is guaranteed to be optimal when the problem concerned satisfies a sufficient condition, i.e., the Principle of Optimality (Bellman, 1953; Puterman, 2014), which states “an optimal policy has the property that whatever the initial state and initial decision are, the remaining

decisions must constitute an optimal policy with regard to the state resulting from the first decision” (Bellman, 1953). Throughout this section, we focus our attention on discrete time systems. For continuous time dynamic systems, the readers are referred to the Hamilton-Jacobi-Bellman equations in optimal control (see, for example, Bertsekas, 2012a).

In particular, for a finite time horizon problem, the decisions are made over a number of stages or decision epochs, denoted by $t = 0, ..., T - 1$. At each decision epoch, after observing the current system state $x_t$ (comprised of one or more information variables that characterise how the system progresses), an action $a_t$ is taken that leads to an immediate reward (cost) of $r_t(x_t, a_t, w_t)$, where $w_t$ is the random disturbance at time $t$ with a known probability distribution. The system then evolves to state $x_{t+1}$ at the next decision epoch, following the transition function $x_{t+1} = f_t(x_t, a_t, w_t)$ with the transition probability $p_t(x_{t+1}|x_t, a_t, w_t)$. After the last decision is made at epoch $T - 1$, the system evolves to $x_T$ in the terminal stage with the salvage value $r_T(x_T)$. The objective of the problem is to find a policy $\pi$, or a sequence of actions ($a_0, a_1, ..., a_{T-1}$) prescribed by $a_t = \pi(x_t)$, that maximises (minimises) the total expected reward (cost) across the entire time horizon. Note that for the expected total reward optimisation criterion (or additive reward functions) the Principle of Optimality is always satisfied (Puterman, 2014). To avoid the technical subtleties, in what follows we focus on discrete state space $S$ and action space $A$, and assume the random disturbance at an epoch is independent of those in the previous epochs. Define the dimension of the state space $S$ as the number of the information variables in the state. The mathematically inclined readers are referred to Puterman (2014) for discussions on more general situations. Before proceeding, it is worth mentioning that when the random disturbance $w_t$ takes only a single value, the problem reduces to a deterministic problem. Perhaps the two most well known deterministic sequential decision-making problems solvable by DP are the Shortest Path problem (Dreyfus, 1969) and the Knapsack problem (Kellerer et al., 2004).

Under the Principle of Optimality, the above-mentioned problem can be solved by backward induction. Denote by $V_t(x_t)$ the value function, or the optimal expected value-to-go from state $x_t$ at epoch $t$ until the end of the time horizon. The value function (for maximisation problems) satisfies the following optimality equations (or Bellman Equations, see e.g., Puterman, 2014),

$$V_t(x_t) = \max_{a_t \in A} \mathbb{E}[r_t(x_t, a_t, w_t) + V_{t+1}(f_t(x_t, a_t, w_t))], \forall x_t \in S, t = 0, ..., T - 1, \quad (1)$$

with the boundary condition $V_T(x_T) = r_T(x_T)$. By recursively solving the optimality equations from the
last stage backwards to time zero, we obtain the optimal value functions and, at the same time, an optimal policy. For this method to work, however, at each stage one has to solve the value function for all states before proceeding to the previous stage. For problems with high dimensional state variables, the solution via this method is simply not practical due to the prohibitive amount of computational time and memory required. The recent development on DP research has been essentially trying to overcome this so called *curse of dimensionality* (Powell, 2011), which is discussed in the last paragraph of this section.

Many sequential decision-making problems in practice do not have a natural termination stage, leading to a rich body of literature studying infinite horizon problems, for which the total expected reward becomes unbounded as the time horizon tends to infinity. To this end, two alternative criteria have been widely used in the literature (Puterman, 2014; Bertsekas, 2012a). The first one applies a discount factor between 0 and 1, say $\beta$, to the future reward, which can be understood as the depreciation of monetary values over time. The total discounted reward is well defined as it is bounded by the sum of a decreasing infinite geometric sequence. In situations where discounting is not appropriate, a meaningful criterion is to consider the long run average reward, or the reward rate per stage. Assuming a stationary system (in which the transition function/probability, the reward function, and random disturbance do not change over time), the Bellman Equations for the total discounted reward criterion take the following form:

$$V(x) = \max_{a \in A} \mathbb{E}[r(x, a, w) + \beta V(f(x, a, w))], \forall x \in S, \quad (2)$$

where the value function $V(x)$ is the optimal discounted value-to-go from state $x$ over an infinite time horizon. Note that there are no more boundary conditions. There is no more dependency on time either under the assumption of stationary systems, which is often satisfied in practice (Bertsekas, 2012a). When such an assumption is not satisfied, a periodic or cyclic DP can be developed (Li et al., 2022a). For brevity we do not include the Bellman equations for the long run average reward criterion but direct the readers to Bertsekas (2012a) and Puterman (2014).

There are mainly three solution algorithms (Tijms, 1994; Puterman, 2014) for infinite horizon problems. The most widely used and understood algorithm is value iteration, or successive approximations as it was called in the early days. Starting from an arbitrary bounded value function vector (e.g., $V_0(x) = 0, \forall x \in S$), this method iteratively updates value functions via the recursive equation below until the successive gaps between iterations

$k+1$ and $k$ are within a predefined threshold.

$$V_{k+1}(x) = \max_{a \in A} \mathbb{E}[r(x, a, w) + \beta V_k(f(x, a, w))], \quad (3)$$
$$\forall x \in S.$$

An alternative algorithm is policy iteration, which starts with an arbitrary policy and then iteratively improves it until no further improvements are possible. Each iteration includes two steps: firstly the expected value-to-go under the current policy is evaluated via a system of equations similar to (2) but for the actions prescribed by the policy; after that a policy improvement step is undertaken to find an improved action for each state that leads to a better value-to-go (Puterman, 2014). In the last algorithm, the system of Bellman Equation (2) are reformulated into a vary large scale linear program, which has one decision variable for each state and one constraint for each state-action pair. Regardless of the solution algorithms, just as in finite horizon problems, the curse of dimensionality remains the biggest hurdle for the implementation of DP.

Various approximation methods have been proposed to improve the scalability of DP, leading to an important and thriving research field called Approximate Dynamic Programming (ADP). According to Bertsekas (2012a), most of the ADP approaches fall into either the value space or policy space. We concentrate on the approaches in the value space (see also §2.21) while we direct readers to Bertsekas (2012a) for the policy space counterparts. The basic idea of the value space approaches is to develop efficient methods to approximate the value functions or the expected value-to-go for a given policy. The most studied methods approximate the value functions via a linear or nonlinear combination of a set of handcrafted feature vectors (functions of the state) weighted by a set of parameters, which are calibrated by a suitable method (Bertsekas, 2012b; Ding et al., 2008). Feature vectors are not always available, in which case Neural Networks have been used to construct feature vectors automatically (Powell, 2011; Bertsekas, 2012a; He et al., 2018). Decomposition is also a popular method, which decomposes the original problem into a number of sub-problems each of which has a much smaller state space and can be solved efficiently by the exact algorithms mentioned above. The assembly of the value functions of these sub-problems provides an approximation to the original value functions (Kunnumkal & Topaloglu, 2010; Li & Pang, 2017). A distinct decomposition approach is Whittle's Restless Bandit framework (Whittle, 1988; Glazebrook et al., 2014; Li et al., 2020), which decomposes the original problem via Lagrangian relaxation, calculates a state dependent index value for each sub-problem and uses these index values directly to derive policies for the original
<figure><img src="image_34.png" /></figure>

problem. Another method in the value space approximates the value functions of a specific policy via Monte Carlo simulation (Chang et al., 2007; Bertsekas, 2012b), which are then used to find an improved policy. An alternative method is called *Q-Learning* (Sutton & Barto, 2018), which approximates the *Q-factor* for each state-action pair. The *Q-factor* for (x, a) is the expected value-to-go by taking action a at state x and then following either a given policy or the optimal policy thereafter. Due to the large number of combinations of state-action pairs, *Q-Learning* is more suitable for problems with a small state space (Bertsekas, 2012b). For an in-depth account on ADP we refer to two seminal books of Powell (2011) and Bertsekas (2012b).

## 2.10. Forecasting¹⁶

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 & Jenkins, 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 & 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
<figure><img src="image_35.png" /></figure>

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 coherence (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 & 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 & Benítez, 2012). A more recent approach to forecast selection is based on the concept of representativeness (Petropoulos & 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 & 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 & Wright, 1999), structured analogies (Green & Armstrong, 2007), prediction markets (Wolfers & Zitzewitz, 2004), and interaction groups (Van de Ven & 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 & 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
<figure><img src="image_36.png" /></figure>

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 & 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 https://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 theory¹⁷

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 & Sigmund, 1998; Cressman, 2003; Sandholm, 2010) and mean-field games, (Huang et al., 2003, 2006, 2007; Lasry & Lions, 2006a, 2006b, 2007; Gomes & 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
<figure><img src="image_37.png" /></figure>

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 & Neyman, 1981; Jaśkiewicz & Nowak, 2018a, 2018b). In a multistage game, the players share the control of a discrete-time dynamic system (state equations) observed over stages (Başar & Olsder, 1999; Engwerda, 2005; Krawczyk & 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 architecture 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 = \{1, ..., n\}$ 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(s_1, ..., s_n) : \prod_{i \in I} S_i \to \mathbb{R}$, where $s = (s_1, ..., s_n)$. Assuming the players are maximisers, the strategy profile $s^N = (s_1^N, ..., s_n^N)$ is a Nash equilibrium if $g_i(s_1^N, ..., s_n^N) \ge g_i(s_1^N, ..., s_{i-1}^N, s_i, s_{i+1}^N, ..., s_n^N)$ 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 = \underset{s_i \in S_i}{\operatorname{argmax}} g_i(s_1^N, ..., s_{i-1}^N, s_i, s_{i+1}^N, ..., s_n^N)$. 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
<figure><img src="image_38.png" /></figure>

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 and Leitmann (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 & Petrosyan, 2018; Petrosyan & Zaccour, 2018).

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

## 2.12. Graphs and networks¹⁸

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 $O(m + n \log n)$ time. For graphs with possible negative arc costs, in $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 & Fulkerson, 1962). The method of Ford and Fulkerson (1957) is pseudopolynomial 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 & 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 & Tarjan, 1988). A *pseudoflow* algorithm for the
<figure><img src="image_39.png" /></figure>

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 $P$, however most of the network optimisation problems that are relevant for practical applications are $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 $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 & Rao, 1994). Exact solution methods for the STP are based on a sophisticated combination of: *reduction techniques* (Gamrath et al., 2017; Rehfeldt & 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 & 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 achieved¹⁹. Facet-defining inequalities are investigated in Padberg and
<figure><img src="image_40.png" /></figure>

Rinaldi (1990b); Jünger et al. (1995). MIP formula-
tions, including the famous *subtour-elimination con-
straints* 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 approxima-
tion algorithms for the TSP see Traub (2020).
Helsgaun (2000)²⁰ provides an efficient implementa-
tion 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
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), instan-
ces 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 vertices²².

**2.13. Heuristics<sup>23</sup>**

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 meth-
ods, by using different techniques ranging from sim-
ple 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 oper-
ational 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 lit-
erature guide us for the methods employed in heu-
ristics. 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 induc-
tion, 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 & 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 & Martí 2013). One way of
doing this is to identify common attributes of the
optimal solution and search among the solutions hav-
ing these attributes only (Glover, 1977). Another

restriction can be applied by eliminating infeasible
solutions considering a combination of decision vari-
ables which dictates incompatible values. Beam
Search (Morton & Pentico, 1993) is a good example
of this class of heuristics which works with a trun-
cated 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 & 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 pro-
posed by Hax and Meal (1973) considers the organ-
isational 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 deci-
sions 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 decom-
position 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 mathemat-
ical models and utilise different strategies to make
the problem tractable which results in a reduced
<figure><img src="image_41.png" /></figure>

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 & 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 & 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 from 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 & 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 & 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
<figure><img src="image_42.png" /></figure>

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. Solyali 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 programming²⁴

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 \ge 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 & Steiglitz, 1982). Since
<figure><img src="image_43.png" /></figure>

$\min(cx) = -\max(-cx)$, one can minimise a "cost" as well maximise a "profit". Linear programs come in pairs: $\{\text{minimise } cx, \text{ subject to } Ax = b, x \ge 0\}$ and $\{\text{maximise } yb, \text{ subject to } yA \le 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 \le 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$ nonzero 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 $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 & 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 \ge 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
<figure><img src="image_44.png" /></figure>

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*⁻¹, and this is repeated with a new scaling bijection (see Strang, 1987; Goldfarb & Todd, 1989; Fang & Puthenpura, 1993; Winston & 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 $\mathcal{P}$, ILP is not, unless $\mathcal{P} = \mathcal{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, 1965b; Grötschel et al., 1981, 1988; Cook et al., 1998). Furthermore, tackling $\mathcal{NP}$-complete problems has benefited greatly from this approach (see, e.g., Dantzig et al., 1954; Applegate et al., 2007).

## 2.15. Integer programming²⁵

Mixed-integer programming (MIP) is an $\mathcal{NP}$-hard generalisation of linear programming (LP; $\S2.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 & 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
<figure><img src="image_45.png" /></figure>

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 & 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 & Wolsey, 1986), which eventually led to effective cut-and-branch algorithms for such problems (e.g., Van Roy & 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 & Schrijver, 1991; Sherali & 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 *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 witnessed 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 & 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 & 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 & Sağlam, 2013).

## 2.16. Nonlinear programming<sup>26</sup>

Nonlinear programming is a generalisation of linear programming (§2.14), in which the objective function or the constraints can be given by general