Data driven optimisation models for the operations

planning in industrial supply chains

Benoit Loger

To cite this version:

Benoit Loger. Data driven optimisation models for the operations planning in industrial supply chains.
Operations Research [math.OC]. Ecole nationale supérieure Mines-Télécom Atlantique, 2023. English.
�NNT : 2023IMTA0389�. �tel-04497849�

HAL Id: tel-04497849

https://theses.hal.science/tel-04497849
Submitted on 11 Mar 2024

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
T HÈSE DE DOCTORAT DE

L’ÉCOLE NATIONALE SUPÉRIEURE

MINES-TÉLÉCOM ATLANTIQUE BRETAGNE
PAYS-DE-LA-LOIRE - IMT ATLANTIQUE

É COLE D OCTORALE N O 648

Sciences pour l’Ingénieur et le Numérique
Spécialité : Informatique

Par
Benoit LOGER
Modèles d’optimisation basés sur les données pour l’optimisation
des opérations dans les Supply Chain industrielles

Thèse présentée et soutenue à IMT Atlantique, Nantes, le 21/12/2023

Unité de recherche : Laboratoire des Sciences du Numérique de Nantes (LS2N)
Thèse N° 2023IMTA0389

Rapporteurs avant soutenance :

André ROSSI Professeur à Université Paris Dauphine
Céline GICQUEL Maître de Conférences à Université Paris Saclay

Composition du Jury :
Président : Nadjib BRAHIMI Professeur à Rennes School of Business
Examinateur : André ROSSI Professeur à Université Paris Dauphine
Céline GICQUEL Maître de Conférences à Université Paris Saclay
Stefan MINNER Professeur à Technical University of Munich
Evren SAHIN Professeure à CentraleSupélec
Dir. de thèse : Fabien LEHUÉDÉ Professeur à IMT Atlantique
Co-encadrants : Guillaume MASSONNET Maître Assistant à IMT Atlantique
Alexandre DOLGUI Professeur à IMT Atlantique
 R ÉSUMÉ EN FRANÇAIS

Dans notre société, la majorité des biens de consommation est produite et distribuée
par des entreprises dont l’objectif est de satisfaire efficacement (à moindre coût et avec une
haute qualité de service) la demande de leurs clients. Depuis de nombreuses décennies, les
logiques de compétitivité ont poussé de nombreuses entreprises à se spécialiser dans une
unique tâche (de production ou de service) et à collaborer avec d’autres entités spécial-
isées. De ce fait, les circuits de production et de distribution sont devenus de plus en plus
complexes, incluant une grande diversité d’acteurs amenés à collaborer pour satisfaire les
besoins finaux. De tels réseaux d’entreprises sont appelés chaînes d’approvisionnement (ou
Supply Chain en anglais). Leur diversification et leur complexification durant les dernières
décennies a poussé les gestionnaires de ces réseaux à caractériser formellement leurs dy-
namiques afin d’en optimiser le fonctionnement. L’identification de ces besoins et prob-
lématiques, génériques commes spécifiques, ont motivé de nombreux travaux de recherches
et de nombreux livres portant sur le management des chaînes d’approvisionnement (Stadtler
et al. (2015),Christopher (2005),Dolgui and Proth (2010),Simchi-Levi (2008)). Ces travaux
ont ainsi contributé à décrire la structure et les dynamiques de ces chaînes d’approvision-
nement en détail, à identifier leurs failles et les potentiels leviers d’améliorations. Cette
thèse s’inscrit dans la lignée de ces travaux et se concentre sur l’optimisation des processus
au niveau opérationnel (planification de production, gestion de stock, ...).
Comme évoqué plus haut, la collaboration entre les différents maillons d’une chaîne
d’approvisionnement est une des clés clairement identifiées pour améliorer sa qualité de
service et diminuer les inefficiences Dolgui and Proth (2010). Alors que l’accroissement des
flux de données industrielles devraient permettre à chaque entreprise de baser ses décisions
sur une vision globale de la chaîne d’approvisionnement, le partage d’information se heurte
bien souvent aux réticences de groupes parfois en compétition et soucieux de préserver leur
avantage conccurentiel. Ainsi, on constate que de nombreux dirigeants et/ou gestionnaires
ne partagent que partiellement les informations nécessaires aux autres membres d’une
chaîne logistique pour adapter leurs décisions. Ce phénomène est à l’origine d’un intérêt
croissant pour la conception de méthodes d’aide à la décision basées sur une connaissance
partielle de l’état de la chaîne d’approvisionnement, incluant les incertitudes concernant

2
le comportement des autres acteurs et des clients finaux.

Contexte scientifique :
Un des problèmes d’optimisation les plus référencés dans la littérature en gestion de
chaîne d’approvisionnement consiste à planifier simultanément l’approvisionnement en
matières premières ou en composants et les opérations de production. Plus particulière-
ment, le problème de dimensionnement de lot (ou lot sizing) et ses nombreux dérivés a
été l’objet de nombreuses études, dont un nombre non négligeable propose des modèles
qui intègrent des incertitudes. Pendant longtemps, l’essentiel de ces travaux se concen-
trait sur des approches de programmation stochastique visant à optimiser la performance
moyenne des décisions et reposant sur des hypothèses restrictives concernant la distribu-
tion des paramètres incertains. Cependant, dans un contexte applicatif, il est souvent
impossible d’obtenir une caractérisation précise de ces distributions et ces approches peu-
vent conduire à proposer des solutions aux performances très variables, notamment lors
de la réalisation de scenarios peu probables susceptibles d’entraîner des perturbations
extrêmement défavorables.
Dans les deux dernières décennies, de nombreux travaux ont donc cherché à pro-
poser des approches d’optimisation permettant d’obtenir des solutions robustes face aux
aléas, possiblement sous-optimales mais dont la qualité reste satisfaisante même lorsqu’on
s’affranchit des hypothèses restrictives concernant la valeur des paramètres incertains. Ce
paradigme est connu sous le terme d’Optimisation Robuste. Si elle prend racine dans
les travaux de pionniers comme Scarf (1958) ou Soyster (1973), de très nombreux mod-
èles et méthodes de résolution ont été proposés dans la littérature depuis la fin des an-
nées 90 et se démarquent notamment de l’approche par optimisation stochastique par
l’absence d’hypothèse sur la distribution des paramètres incertains inclus dans les mod-
èles d’optimisation. Pour remplacer la description du comportement aléatoire de ses
paramètres, un modèle d’optimisation robuste considère l’ensemble de leurs valeurs pos-
sibles, appelé ensemble d’incertitude. La définition et la formulation de tels ensembles est
au centre de la recherche en optimisation robuste car le choix d’un ensemble influence à
la fois la qualité des solutions obtenues et la capacité de méthodes de résolution classique
telles que la programmation linéaire pour résoudre ces modèles en pratique.
Dans cette thèse, nous nous intéressons plus spécifiquement à des méthodes d’optimisa-
tion robuste dirigées par les données, un domaine de recherche particulièrement actif
au cours des deux dernières décennies. Durant cette période, on assiste en effet à une

3
politique généralisée d’automatisation des opérations de production et de distribution,
accompagnée d’une collecte généralisée des données de mesure relatives aux processus
industriels. L’intégration de cette information aux outils d’aide à la décision, grâce no-
tamment à des approches statistiques ou d’apprentissage automatisé, permet souvent
d’améliorer significativement la qualité des solutions proposées. Dans cette dynamique,
l’optimisation robuste s’avère être une approche particulièrement bien adaptée, à même
d’intégrer directement ces données dans la caractérisation des ensembles d’incertitude
considérés. Cependant, une étude de la littérature montre que l’application de ces ap-
proches aux problèmes propres à la gestion de la chaîne d’approvisionnement n’a été que
très peu étudiée. En effet, la plupart de ces ensembles d’incertitude mènent à des modèles
d’optimisation de grande taille dont la résolution par des approches classiques de program-
mation linéaire est souvent trop coûteuse en temps et en mémoire, limitant leur champs
d’application à des problèmes simples n’intégrant que peu de paramètres incertains.
Les travaux de cette thèse ont donc pour objectif de proposer des approches d’optimisa-
tion robuste dirigées par les données pour des problèmes d’approvisionnement et de pro-
duction. Nous nous intéressons tout particulièrement à la qualité des solutions proposées
par nos modèles en comparaison des modèles déterministes (sans incertitudes), de modèles
robustes plus classiques et de modèles d’optimisation stochastique. Nos contributions por-
tent notamment sur la formulation des ensembles d’incertitude induits et sur la résolution
des modèles qui en découlent.

Structure du manuscrit :
Le manuscrit est composé des 6 chapitres suivants :
Le Chapitre 1 introduit le concept d’optimisation robuste ainsi que les modèles et les
méthodes de résolution les plus utilisés dans la littérature. De nombreux travaux an-
térieurs à cette thèse ont proposé des ensembles d’incertitude permettant de modéliser
plus précisément les incertitudes afin de les intégrer aux modèles d’optimisation robuste.
Une revue de la littérature est donc faite pour présenter ces différentes approches et
notamment les approches dirigées par les données. Pour introduire plus précisément le
contexte applicatif de cette thèse, ce chapitre présente également une revue de littérature
des applications de l’optimisation robuste pour des problèmes de gestion de la chaîne
d’approvisionnement.
Le Chapitre 2 présente les différentes contributions de cette thèse et les positionne vis-
à-vis de l’état de l’art.

4
Le Chapitre 3 présente une méthode d’approximation pour une classe d’ensembles
d’incertitude dirigés par les données. Ces ensembles sont obtenus à l’aide d’un algo-
rithme d’apprentissage (connu sous l’appellation Support Vector Clustering) et les utiliser
dans des modèles d’optimisation robuste permet d’obtenir des solutions de meilleure qual-
ité que les modèles classiques. Cependant, leur formulation complexe augmente de façon
significative le temps nécessaire à l’obtention de solutions optimales et limite, de fait, leur
champs d’application à des problèmes de petite taille. L’approximation introduite dans ce
chapitre réduit ce temps de calcul de façon significative et permet d’obtenir des solutions
optimales aux problèmes d’optimisation robuste en un temps raisonnable.
Le Chapitre 4 est consacré à un cas d’application pratique inspiré du domaine aéronau-
tique. Le problème consiste à planifier les opérations d’approvisionnement et de produc-
tion d’une ligne d’assemblage dans le but de satisfaire les besoins en différents produits
finaux. Les opérations d’approvisionnement sont soumises à des incertitudes émanant de
la collaboration avec un prestataire logistique et menant à une diminution de l’efficacité
du système de production. Ce chapitre introduit deux modèles d’optimisation robuste
permettant de garantir un haut niveau de satisfaction de la demande en produits finis.
Le premier modèle utilise une approche classique de la littérature et le second utilise
l’algorithme de support vector clustering présenté au chapitre précédent, dans le but
d’obtenir des solutions optimales en un temps raisonnable. Les expérimentations mon-
trent l’intérêt des modèles robustes en comparaison à une approche déterministe, ainsi
que la supériorité des solutions obtenues avec le modèle dirigé par les données.
Le Chapitre 5 est consacré à un problème plus général d’approvisionnement face à des
incertitudes sur la demande et sur le délai de livraison. Ce chapitre introduit un nou-
vel ensemble d’incertitude dirigé par les données qui, comparé à un modèle robuste plus
classique, permet d’obtenir des solutions de meilleure qualité dans les cas de corrélation
étudiés. Cet ensemble d’incertitude diffère des modèles introduits jusqu’alors. En ef-
fet, il est basé sur l’hypothèse que la demande suit un processus stochastique bien défini
(mais dont les paramètres sont incertains). Nous utilisons ici un modèle ARMA (Au-
toregressive Moving Average), communément utilisé dans des méthodes prédictives, pour
définir les bornes de notre ensemble d’incertitude. Les expérimentations montrent que
cette approche permet d’obtenir des solutions moins sensibles aux scénarios de demande
extrêmes en comparaison d’une approche d’optimisation stochastique. De plus, comparé
à une approche d’optimisation robuste classique de la littérature, cette approche permet
d’obtenir de meilleures performances en moyenne et n’induit pas d’augmentation du temps

5
de résolution des modèles.
Le Chapitre 6 résume les conclusions de cette thèse et introduit diverses perspectives
pour de futurs travaux de recherche.
Les travaux présentés dans cette thèse ont mené à la soumission d’un article scientifique
en cours de publication dans une revue internationale ainsi qu’à des présentations dans
diverses conférences nationales et internationales. Deux autres articles sont également en
cours de rédaction.

6
 ACKNOWLEDGEMENT

Je souhaite remercier chaleureusement mes directeurs Fabien Lehuédé, Guillaume Mas-

sonnet et Alexandre Dolgui pour leur soutien tout au long de ces trois années. Je tiens
à leur exprimer ma plus profonde gratitude pour la confiance qu’ils m’ont accordée, pour
leurs conseils, leur écoute, et plus généralement, pour tous les échanges que nous avons
pu avoir (scientifique ou non).

Je remercie Céline Gicquel et André Rossi d’avoir accepté la charge d’être relecteur
de cette thèse et pour leurs remarques qui ont abouti à cette version finale du manuscrit.
Je remercie également Nadjib Brahimi, Evren Sahin ainsi que Stefan Minner pour leur
participation en tant que membre du jury et pour nos échanges qui offrent de nouvelles
perspectives pour de futurs travaux. Je remercie tout particulièrement Stefan Minner
pour son accueil au sein de l’université technique de Munich et pour nos échanges lors de
notre collaboration.

Je remercie également l’ensemble de l’équipe Modelis d’IMT Atlantique de m’avoir

offert un cadre de travail convivial tout au long de ma thèse. En partageant avec moi un
bureau, un café, une bière ou une cible d’aimant, vous avez su me remotiver lorsque je
rencontrais des moments difficiles.

La rédaction de cette thèse n’aurait pas été possible sans le soutien de mes parents,
Dominique et Gilles, de mes frères et sœurs, Guillaume et Hélène, et de mes plus proches
amis, ils se reconnaîtront. Enfin, je remercie ma copine Mélissa, sans qui rien n’aurait été
possible. Elle a été mon plus grand soutien avant et pendant la rédaction de cette thèse,
ma plus grande force et ma plus grande fierté.

7
 TABLE OF C ONTENTS

Introduction 13
1 Scientific context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Overview of the thesis and notations . . . . . . . . . . . . . . . . . . . . . 15

1 Literature Review 17
1 Robust Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.1 Solving RO models . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.1.1 Dual reformulation with polyhedral uncertainty set . . . . 20
1.1.2 Iterative cutting plane methods . . . . . . . . . . . . . . . 22
1.2 Constructing Uncertainty sets . . . . . . . . . . . . . . . . . . . . . 23
1.3 Data-Driven uncertainty sets . . . . . . . . . . . . . . . . . . . . . 25
1.3.1 Gaussian Mixture Models . . . . . . . . . . . . . . . . . . 27
2 Robust inventory management . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1 Lot-sizing models . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Static robust inventory and production management . . . . . . . . 31
2.3 Dynamic robust inventory and production management . . . . . . 33

2 Contributions 37

3 Approximate kernel learning uncertainty set 41

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2 Support vector clustering-based uncertainty set . . . . . . . . . . . . . . . 43
3 Approximation of the SVC-based uncertainty set . . . . . . . . . . . . . . 46
3.1 Initial support vectors selection by K-medoid clustering . . . . . . . 47
3.2 Procedure OPT -α . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3 Reduction of the number of support vectors . . . . . . . . . . . . . 49
3.4 Summary of the approach . . . . . . . . . . . . . . . . . . . . . . . 50
4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.1 Uncertainty set comparison . . . . . . . . . . . . . . . . . . . . . . 52
4.2 Application to optimization problems . . . . . . . . . . . . . . . . . 54

9
 4.2.1 Robust 0-1 Knapsack Problem . . . . . . . . . . . . . . . 55
4.2.2 Robust Lot-sizing Problem . . . . . . . . . . . . . . . . . 57
4.2.3 Robust Traveling Salesman Problem . . . . . . . . . . . . 61
5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4 A Data-Driven Robust Optimization approach to supply planning 65

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
2 Problem statement and deterministic formulation . . . . . . . . . . . . . . 67
3 Robust models for uncertain setup times . . . . . . . . . . . . . . . . . . . 69
3.1 Budget-based robust formulation . . . . . . . . . . . . . . . . . . . 69
3.2 SVC-based robust formulation . . . . . . . . . . . . . . . . . . . . . 71
4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.1 Scalability of the SVC-based model . . . . . . . . . . . . . . . . . . 74
4.2 Comparing U Γ and U AC . . . . . . . . . . . . . . . . . . . . . . . . 76
5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5 Data-Driven Robust inventory management under uncertain demand

and lead time 81
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
2 Problem and robust solution approach . . . . . . . . . . . . . . . . . . . . 83
2.1 Deterministic model . . . . . . . . . . . . . . . . . . . . . . . . . . 84
2.2 Robust optimization model . . . . . . . . . . . . . . . . . . . . . . 85
3 Data-driven uncertainty set . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.1 Data driven uncertainty set based on ARMA(p,q) process . . . . . 87
3.1.1 General ARM A(p, q)-based uncertainty set . . . . . . . . 87
3.1.2 Examples for special cases of demand processes . . . . . . 88
3.2 Joint demand and lead time uncertainty set . . . . . . . . . . . . . 89
3.3 A constraint generation approach . . . . . . . . . . . . . . . . . . . 91
4 Computational experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.1 Sample average approach . . . . . . . . . . . . . . . . . . . . . . . . 93
4.2 Uncertain correlated demand . . . . . . . . . . . . . . . . . . . . . 94
4.3 Uncertain demand-dependent lead times . . . . . . . . . . . . . . . 98
5 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6 Conclusions and Perspectives 103

10
Conclusion 103

Appendices 109
A Support Vector Clustering and Kernel function . . . . . . . . . . . . . . . 109

Bibliography 111

11
 I NTRODUCTION

In today’s society, the production and distribution of goods is largely carried out by
large companies, which are responsible for satisfying the demand or needs of end cus-
tomers. In the current context, where international trade is encouraged and the economic
competitiveness of companies is a priority, more and more companies are specializing in a
single activity, collaborating with other specialized companies or subcontracting some of
the tasks they used to carry out. As a result, production and distribution channels have
become increasingly complex, incorporating more and more organizations whose collab-
oration is necessary to ensure that market needs are met. Such networks of companies
are often referred to as Supply Chains (SC). Their diversification and expansion over the
past decades has motivated numerous research studies and books focused on Supply Chain
Management (SCM) (Stadtler et al. (2015); Christopher (2005); Dolgui and Proth (2010);
Simchi-Levi (2008)), with the objective to describe them in detail, identify their bottle-
necks and opportunities of improvement. By searching through the literature on SCM, we
can easily see that it regroups a wide variety of problems going from the more strategic
(e.g. defining the location of production sites and warehouses, selecting suppliers, etc.)
to the more operational considerations (e.g. scheduling of daily production, transmission
of orders to a supplier, etc.). This thesis focuses mainly on the operational level, we refer
the interested reader to the wide literature on SCM for more detail on the strategic part.
Many definition of what is a supply chain co-exist because of the wide range of organi-
zations that this term aims to describe. Christopher (2005) defines a supply chain as a « ...
network of organizations that are involved, through upstream and downstream linkages, in
the different processes and activities that produce value in the form of products or services
in the hand of the ultimate customer. ». In other words, a SC is composed of several
organizations who collaborate through material, financial and information exchange.
For instance, think of a production firm that collaborates with a Logistic Service
Provider (LSP) in charge of the reception and storage of raw materials and their delivery
to various locations on the production line. In this simple example, there is an upstream
material flow (i.e. the LSP deliver raw material on the production line) and downstream
financial and information flows. The financial flow is the price paid by the production

13
firm for the services provided and the information flow is the transmission of the quantity
of raw materials needed on the production line.
One point of agreement over all definitions of SC is that its main objective is the satis-
faction of the final customer demand. To achieve this goal, ensuring a good collaboration
of all actors of a SC has been identified as a key ingredient to ensure the reliability of
material flows. The quality of this collaboration may rely on information sharing, giving
access to the state of the SC as a whole to each organisation in order to improve their
decision making with respect to a global target, rather than limiting their knowledge
to partial information from their direct upstream and downstream collaborators (Barlas
and Gunduz (2011)). However, complete collaboration is hard to implement in practice
because each company involved in a SC also has to maintain its competitiveness. For
example, a manufacturer selling its products to several retailers may share the level of
congestion on its production line to inform them about the expected lead time of their
orders. Such additional knowledge may in turn lead the retailers to adapt their decisions
to reduce the risk of shortage, for instance by anticipating their needs or requesting a spe-
cial delivery from another producer. While this common knowledge is likely to improve
the service level of the whole SC, the manufacturer may perceive it as an additional risk
of losing contracts, making it harder to prioritize orders coming from a particular retailer.
Such issues have lead SC managers to only partially disclose the information on their own
activity to other members, raising the interest for decision support tools that are able
to base their optimization procedures on a partial knowledge about the state of the SC,
including uncertainties on the behavior of other actors or the end customers.

1 Scientific context
In most applications, SCM involves solving production planning problems consider-
ing simultaneously the production, procurement and distribution decisions (Pochet and
Wolsey (2006)). Among the classical SCM research topics, Lot Sizing Problems (LSP)
have been the focal point of a large part of the literature, counting numerous contri-
butions since the pioneering works of Harris (1913) and Wagner and Whitin (1958). A
rising challenge that recent works have attempted to address is related to problems in
which uncertainty occurs, with the objective to mitigate the impact of uncertainties on
the reliability of production plan and/or production costs.
In this context, Stochastic Programming has been one of the most popular approaches

14
and studies have considered a wide variety of problems and sources of uncertainties (see
e.g. Brahimi et al. (2017); Metzker Soares (2022)). In order to integrate uncertainties
in the LSP, stochastic programming approaches often rely on strong assumptions on the
distributions of uncertain parameters. In addition, the objective of such models is to
minimize the expected cost of solutions thus taking the risk to face extreme costs or to
produce unreliable production plans in extreme cases of uncertain parameters. Thus,
in order to reduce the risk of bad performances in extreme conditions, many decision
makers seek for a "robust" solution (i.e. a solution that performs reasonably well for all
realizations of uncertain parameters).
While the last two decades have witnessed the rapid development of Robust Optimiza-
tion (RO) models for various problems, its introduction for SCM applications dates back
to the work of Scarf (1958); Gallego (1992); Moon and Gallego (1994). RO models have
the advantage that the approach does not require specific assumptions on the distribution
of uncertain parameters, since they only consider that some of the problem parameters
can take any value within a bounded set, often called uncertainty set. A solution is then
said to be robust if it remains feasible for all values of uncertain parameters within this
set.
In parallel, the volume of data collected by the actors of supply chains grows contin-
uously, giving access to large databases of historical values of these uncertain parameters
(e.g. demand, lead time, picking time) to improve the decision making. It turns out that
many RO methods offer a suitable framework to incorporate directly these historical data
in the formulation of optimization models (Bertsimas et al., 2017), which has led to a new
research stream referred to as Data Driven Robust Optimization (DDRO).

2 Overview of the thesis and notations

In this thesis, we study problems related to SCM in an uncertain context through the
lens of data-driven robust optimization.Chapter 1 introduces the concept of robust opti-
mization, giving the most common models and solution methods. It presents a structured
literature review on the different approaches to model the uncertainty, going from models
using basic statistical information to data-driven models built with machine learning algo-
rithms. The second part of our survey is dedicated to the application of RO in production
planning and inventory management. It gives an overview of the different problems that
have been studied under the RO framework, the most common models and optimization

15
algorithms used in this area. In Chapter 2, we summarize the contributions presented in
detail in Chapter 3, 4 and 5 of this thesis. Finally, we conclude this manuscript by a sum-
mary of the contributions developed in these chapters, giving some managerial insights
and describing promising perspectives for future research.
Notations: In the rest of the manuscript we use the following notations in order
to simplify mathematical expressions. All scalar parameters or variables are denoted by
lower case letters such as aij . Bold lower case letters such as x denote vectors of variables
or parameters and bold uppercase letters are used to denote matrices such as A whose
j − th row is denoted by aj . We denote by [m] the set of integer values going from 1 to
m and ∥·∥p denote the lp -norm for all p ∈ N+ ∪ {∞}. Sets are denoted by calligraphic
upper case letters such as D, S or B.

16
 Chapter 1

L ITERATURE REVIEW

In this chapter, we introduce the concept of Robust Optimization and review the exist-
ing literature related to uncertainty set construction, Data Driven Robust Optimization
and Robust Supply Chain Management.

1 Robust Optimization
Since the end of 1990’s, Robust Optimization (RO) has become a popular approach
to solve optimization problems involving uncertain parameters. The key aspect of RO
is that it does not rely on distributional assumptions of uncertain parameters unlike
other standard approaches such as Stochastic Optimization (SO) or Stochastic Dynamic
Programming (SDP). The central concept of RO is to consider that some of the problem
parameters aj can take any value within a bounded set Uj , often called uncertainty set,
instead of assuming a fixed one. A solution is then said to be robust if it remains feasible
for all possible values aj ∈ Uj . In this section, we consider the general form of robust
linear optimization problems given by:

max cT x
s.t. aTj x ≤ bj ∀j ∈ [m], ∀aj ∈ Uj (1.1)
x ∈ Rn+

where x is an n-dimensional decision variable vector, c a vector of cost parameters,

aj ∈ Rn are vectors of uncertain parameters of the problem for all j ∈ [m].
The concept of (RO) was first introduced by Soyster (1973) who defines a framework
in which each uncertain parameter aij belongs to the interval [āij − âij , āij + âij ] where
āij is called the nominal value and âij represents the maximum deviation of aij from āij .

17
 aij −āij
Introducing additional notations zij = âij
. The resulting uncertainty set is defined as:

UjBox = {aij = āij + zij âij | ∥zj ∥∞ ≤ 1} (1.2)

. From the definitions (1.1) and (1.2), it follows that any solution using Soyster’s assump-
tions is protected against all possible values of the uncertain parameters within U Box .
While this property may be desirable in some very specific settings, it turns out that
when parameters aij are in fact random variables with support [āij − âij , āij + âij ], this
approach often becomes overly conservative. It is well established that in most applica-
tions, it is highly unlikely that all parameters simultaneously take their extreme values,
and the probability that this situation occurs decreases as the number of random param-
eters increases. In such cases, the solutions obtained using Soyster’s modeling framework
guarantee their robustness against unlikely realizations of the problem parameters, of-
ten at the expense of a significant performance reduction. Such solutions are called
over-conservative. This issue has been a major topic for researchers over the last three
decades, pursuing the right trade-off between the average performance and the protection
level against bad parameters realizations.
The first attempt to reduce over-conservatism was proposed simultaneously by Ben-
Tal and Nemirovski (1998) and El Ghaoui et al. (1998) who studied robust convex and
semidefinite programs and their tractability. They start from the fact that, using the l2 -
norm, the cumulated deviation of uncertain parameters from their nominal values evolve
as an ellipse in the parameter space. Thus, they introduce an ellipsoidal uncertainty
set whose radius is controlled by a parameter Ωj . The use of this parameter provide
an natural interpretation to control the conservatism of RO models as increasing the
value of Ωj would be used to ensure a larger constraint satisfaction level but will often
result in a reduction of performance. Under such uncertainty set, they propose to solve
the robust problem (1.1) by replacing it by an equivalent and tractable deterministic
reformulation (the Robust Counterpart (RC)). Among other results, they prove that the
robust counterpart is a Second Order Cone Program (SOCP) when the initial problem
is a Linear Program (LP). Later on, Ben-Tal and Nemirovski (2000) extended this work
and provide an upper bound on the constraint violation level for robust LP under such
uncertainty sets when parameters are independent and symetrically distributed random
variables. The ellipsoidal uncertainty set is given by:
n o
UjΩ = aij = āij + zij âij ∥zj ∥2 ≤ Ω2j (1.3)

18
 Following their works, Bertsimas and Sim (2003, 2004) introduce a polyhedral uncer-
tainty set whose size is also controlled by a parameter Γj ∈ {0, . . . , n} called budget of
uncertainty. This uncertainty set, often called budget-based uncertainty set is given by
the following expression:

UjΓ = {aij = āij + zij âij | ∥zj ∥1 ≤ Γj } (1.4)

As in the ellipsoidal case, Γj provides a natural natural interpretation to control the

robustness of robust models as Γj represents the maximum number of parameters which
can deviate from their nominal value āij in the i-th constraint. In other words, if Γj = 0,
all uncertain parameters will take their nominal value. On the contrary, if Γj = n they
will all take their worst value. However in this case, the use of ∥·∥1 ensures that the
robust counterpart of a LP remains linear, thus providing a significant computational
advantage compared to the ellipsoidal uncertainty set (1.3). With the help of numerical
experiments on a binary knapsack problem and a linear portfolio management problem,
they prove that their approach reduces the conservatism of robust solutions (compared to
U Box ) while allowing large Mixed Integer Linear Problems (MILP) to be solved efficiently.

As an example, Table 1.1 show the performance reduction and the probability of
constraint violation obtained for the binary knapsack problem under uncertain weights of
Bertsimas and Sim (2003, 2004). We compare the solutions obtained by the deterministic
model and the robust models using U Box and U Γ and we evaluate the empirical probability
of constraint violation on 10000 random scenarios. On this example, we observe that
using Soyster’s method provide 100% of constraint satisfaction but reduce the objective
value by 5.19% in comparison with the deterministic solution. On the other hand, using
the approach of Bertsimas and Sim (2003, 2004) with different budget provide trade-
off solutions and when Γ = 20.0 the risk of constraint violation is only 0.05% and the
objective value is reduced by only 0.75%.

These three uncertainty sets are at the root of much of the work on Robust Opti-
mization that has followed in the last two decades. Their findings introduces the the
central challenges of robust optimization research: producing high quality solutions while
maintaining computational tractability.

19
Table 1.1 – Performance reduction and risk of constraint violation for the binary knapsack
problem of Bertsimas and Sim (2003, 2004) under box and budget based uncertainty sets.

Method Objective value Performance reduction Constraint violation

Deterministic 8135 0.0% 49.79%

Γ=0 8135 0.0% 49.79%
Γ=5 8123 0.14% 25.28%
Γ = 10 8105 0.35% 6.68%
Γ = 15 8085 0.61% 0.74%
Γ = 20 8074 0.75% 0.05%
U Box 7712 5.19% 0.0%

1.1 Solving RO models

Since RO problems as (1.1) involves an infinite number of constraints when parameters

aij are continuous, using a solver directly to compute an optimal solution is inconceivable.
In this section, we quickly review the most common methods encountered in the literature
to tackle RO models.

1.1.1 Dual reformulation with polyhedral uncertainty set

The most common approach in the literature to solve RO models is to formulate a

deterministic and tractable RC using duality theory. As an example, we describe how to
derive a linear robust counterpart from the robust problem (1.1) when the uncertainty
set is the budget-based uncertainty set (1.4) of Bertsimas and Sim (2003, 2004). Let us
consider the following problem:

max cT x
s.t. max aTj x ≤ bj , ∀j ∈ [m] (1.5)
aj ∈Uj

x ∈ Rn+

This problem is equivalent to (1.1) as any feasible solution x satisfying constraint max aTj x ≤
aj ∈Uj
bj remain feasible for all realization ãj ∈ Uj such that ãTj x ≤ max aTj x.
aj ∈Uj

20
Under the budget-based uncertainty set, this problem can be formulated as:

max cT x
n
( n )
X X
s.t. āij xij + max âij xij zij ≤ bj , ∀j ∈ [m] (1.6)
zj ∈[−1,1]n ∥zj ∥1 ≤Γj
i=1 i=1

x ∈ Rn+

Considering an optimal solution x∗ to this problem, for each j ∈ [m], the inner maxi-
nP o
n ∗
mization problem max
n i=1 â ij x ij zij is equivalent to a linear optimization
zj ∈[−1,1] ∥zj ∥1 ≤Γj
problem. As values in x are all positive or null, we can restrict zj ∈ [0, 1]n and thus
replace the l1 -norm ∥zj ∥1 by a simple sum ni=1 zij . This gives us the following equivalent
P

LP:
n
âij x∗ij zij
X
max
i=1
n
X
s.t. zij ≤ Γj (1.7)
i=1

0 ≤ zij ≤ 1, ∀i ∈ [n].

Associating the dual variables qj and rij , i ∈ [n] to the constraints. We now consider
the dual of problem (1.7)
n
X
min qj Γj + rij
i=1

s.t. qj + rij ≥ âij x∗i ∀i ∈ [n] (1.8)

qj ≥ 0
rij ≥ 0 ∀i ∈ [n]

Since the primal problem (1.7) is feasible and bounded, its dual problem (1.8) is also
feasible and bounded and their optimal values coincide. Therefore, we can replace the
j-th inner maximisation problem of (1.6) by the objective function of (1.8), with the

21
associated constraints. This gives us the (linear) robust counterpart of problem (1.1):

max cT x
n
X n
X
s.t. āij xi + qj Γj + rij ≤ bj ∀j ∈ [m]
i=1 i=1

qj + rij ≥ âij xi ∀i ∈ [n], j ∈ [m] (1.9)

xi ≥ 0 ∀i ∈ [n]
rij ≥ 0 ∀i ∈ [n], j ∈ [m]
qj ≥ 0 ∀i ∈ [m]

which can be solved by linear programming solvers. Note that a similar approach can
be applied to obtain the RC of a robust LP under any polyhedral uncertainty set. For
more details on the formulation of the RC under ellipsoidal uncertainty sets, interested
readers can refer to Ben-Tal and Nemirovski (1998); El Ghaoui et al. (1998); Ben-Tal and
Nemirovski (2000).

1.1.2 Iterative cutting plane methods

As an alternative approach, Mutapcic and Boyd (2009) have proposed to solve RO

models using an iterative cutting plane algorithm. Their method rely on solving a relaxed
RO problem of the form:

max cT x
s.t. aTj x ≤ bj , aj ∈ Ωj , ∀j ∈ [m] (1.10)
x ∈ Rn

where Ωj is a subset of Uj . Given an optimal solution x∗ obtained after k iterations, the

algorithm solves the inner optimization problem to find the worst-case parameter values
(k)
aj . If the j-th constraint is violated, a new constraint is added to the formulation (1.10).
Then, if no more violating realization is found, the current solution x∗ is returned. If
more than one constraint is subject to uncertainties the algorithm can solve one inner
problem per constraint (Bertsimas et al. (2016)). If several constraints contain the same
uncertain parameters, one could formulate a single inner problem corresponding to all the
constraints simultaneously (Bienstock and Ozbay (2008); Thorsen and Yao (2017)).
Such approaches have been proposed to solve robust linear problems, convex prob-

22
lems and quadratic problems (Ben-Tal et al. (2015)). The convergence of such algorithm
have been studied under ellipsoidal uncertainty set (Mínguez and Casero-Alonso (2019))
and under general polyhedral uncertainty set (Bertsimas et al. (2016)). Several papers
compared the performances of such methods with a reformulation approach with ellip-
soidal and budget-based uncertainty sets (Fischetti and Monaci (2012); Bertsimas et al.
(2016)). To date, the latest and more complete comparison (Bertsimas et al. (2016))
claim that there is no significant evidence that one approach outperforms the other. In
fact, the method to be preferred seem to depend on the considered uncertainty set (el-
lipsoidal, budget-based) and on the class of problem (LP, MILP). The systematic use
of reformulation to solve robust problems should therefore be questioned, particularly
when introducing new families of uncertainty sets.Based on the same idea, some papers
also propose dedicated iterative algorithms to solve particular problems (Goerigk et al.
(2015); Bienstock and Ozbay (2008); Montemanni (2006)) or study the use of approximate
solution methods to solve the relaxed robust problem or to identify violating realizations
(Pätzold and Schöbel (2020); Thevenin et al. (2022)).

1.2 Constructing Uncertainty sets

Uncertainties on parameters of optimization problems arise for various reasons. Pa-
rameters can model random events such as the demand of final customers of a SC or
correspond to elements for which an exact measurement is not possible, leading to errors.
However, in most models all types of uncertain parameters are considered as random
variables and this will be the case in the rest of this manuscript. As mentioned at the be-
ginning of this chapter, the construction of uncertainty sets is at the center of researchers
considerations because of its significant role in both the quality of the obtained solutions
and the computational tractability of robust models. The practical efficiency of classical
uncertainty sets such as UjBox , UjΓ and UjΩ have been demonstrated on numerous appli-
cations (Aharon et al. (2009); Gabrel et al. (2014); Sözüer and Thiele (2016)). However,
they systematically define the set of possible values for each parameter as an interval cen-
tred on a nominal value with symmetrical and independent variations. Thus, by ignoring
potential asymmetries or correlations in the distributions of uncertain parameters, using
these uncertainty sets can lead to over-conservative or unreliable decisions in practice.
Moreover, the probability bounds on the constraint satisfaction obtained by using such
uncertainty sets are reliable if and only if uncertain parameters are indepently and symet-
rically distributed. For this reason, numerous papers have attempted to design alternative

23
uncertainty sets that are more suitable to model a wider range of behaviours for uncertain
parameters.
Pachamanova (2002) and Bertsimas et al. (2004) generalize the concept of norm-based
uncertainty set to any arbitrary norm as:

Uj = {aj | ∥M (aj − āj )∥ ≤ ∆j } (1.11)

where ∥·∥ can be any given norm, the budget ∆j bounds the cumulative deviation of
uncertain parameters from their nominal value āj and the matrix M is an invertible
matrix used to integrate statistical information in the uncertainty set. As an example,
considering the l2 -norm, if M is the identity matrix In then Uj is a sphere. If M is
a diagonal matrix whose entries are standard deviations of uncertain parameters then
Uj is the classic ellipsoidal uncertainty set. And if M = Σ−1/2 with Σ the covariance
matrix then Uj is an ellipsoid stretched in the direction of greatest dispersion. This last
matrix is a way to integrate correlation information in the construction of uncertainty
sets. Bertsimas et al. (2004) also prove that the robust counterpart of any problem under
a norm based uncertainty set can be obtained from its dual norm. They derive probability
bounds on constraint violation for the general norm-base which coincides with existing
results from El Ghaoui et al. (1998) and Bertsimas and Sim (2003, 2004).
The concept of general norm based uncertainty set was used later on by Chen et al.
(2007) to propose an asymmetric uncertainty set based on backward and forward deviation
− −
measures. Their approach is to define aj = āj + (â+ +
j − âj ) with âj , âj ∈ Rn+ and to
define an assymetric uncertainty set:
n o
− −
Uj = aj aj = āj + (â+ +
j − âj ), ∥P aj + Qaj ∥ ≤ ∆j (1.12)

where P and Q are diagonal matrices built from forward and backward deviation mea-
sures.
More recently Jalilvand-Nejad et al. (2016) proposed an alternative way to integrate
correlation information in the construction of uncertainty sets. They propose a polyhedral
formulation based on the correlation matrix of uncertain parameters. They use experi-
mental results on a robust job shop scheduling problem with uncertain processing times
to show that their approach outperform Pachamanova (2002) under normal marginal dis-
tributions. Daneshvari and Shafaei (2021) develop an alternative formulation to integrate
correlation information and improve the results on the same set of experiments.

24
 As an other approach to reduce conservatism, Bandi and Bertsimas (2012) propose to
use the central limit theorem to tune the so-called budget of uncertainty. Assuming that
uncertain parameters are independent and identically distributed the corresponding set
is:
n √
( )
X
Uj = aij = āij + zij σ | zij | ≤ Γj n (1.13)
i=1

where σ is the standard deviation of uncertain parameters.

Finally, Poss (2013) developped an other way to deal with over-conservatism. He uses
the budget-based approach of Bertsimas and Sim (2003, 2004) in which the budget Γj is
a function of the decision variables and defines an uncertainty set

UjΓ = {aij = āij + zij âij | ∥zj ∥1 ≤ γj (x)} (1.14)

He derives a linear robust counterpart when the function γj (x) is a linear or piecewise lin-
ear increasing function and proves that the above uncertainty set reduces the conservatism
in problem with few non-zero variables. His approach can be applied for problems with
binary and bounded integer or continuous variables. He illustrates the performance of his
approach over the classic budget-based approach on the Robust 0-1 Knapsack problem
and its fractional version.
The variations of uncertainty sets presented above offer promising ways to reduce the
conservatism of robust models in particular cases under specific assumptions. However,
they either rely on strong assumptions on the distributions of uncertain parameters or
necessitate the tuning of parameters (nominal values, maximum deviations, budget pa-
rameters) that strongly influences the efficiency of the resulting robust models. Making
robust optimization a tool of interest for practitioners thus requires uncertainty sets that
are both as flexible as the asymmetric ones of Chen et al. (2007) and parameterized with
simple and intuitive parameters such as the budget of uncertainty introduced in Bertsimas
and Sim (2003, 2004).

1.3 Data-Driven uncertainty sets

Up to now, we illustrated the fact that in order to reduce conservatism, the traditional
approach to include uncertain parameters in robust optimization models relied on offline
procedures and/or assumptions to determine their range of admissible values or distri-
butional characteristics. In this classical framework, upfront data analysis is sometimes

25
employed to define some input parameters to the problem under study. However, with
the emergence of large data sets collected in a variety of applications, the last two decades
have led researchers and decision makers to challenge this approach. A new trend has
emerged, which allows to directly incorporate the data into a mathematical formulation,
often through a combination of statistical analysis and/or dedicated algorithms, and tailor
the extracted knowledge to the problem under study. Within these developments, RO has
emerged as a suitable approach that enables the practitioners to both integrate raw infor-
mation into optimization models while retaining good generalization properties, see e.g.,
Bertsimas et al. (2017). Such models are referred to as Data-Driven Robust Optimization
(DDRO) models in the literature.
In the research literature on DDRO, most efforts focus on leveraging various statistical
and machine learning approaches on historical data to build uncertainty sets that embed
structural properties of the model parameters. One of the major goal of DDRO is to
combine precision in the description of the uncertain parameters and simplicity of use for
practitioners. In this section we review the main existing approaches from the literature
to obtain exploitable data-driven uncertainty sets.

1.3.0.1 Statistical Hypotestis Test In their well-known paper Bertsimas et al.

(2017) use statistical hypothesis-test (SHT) to build uncertainty sets with different as-
sumptions on the distributions (i.e. known finite support, continuous with independent
marginals, none). Each of their uncertainty sets corresponds to a single tractable robust
counterpart formulation (LP, SOCP) and is related to a given bound on the probability of
constraint violation (valid under independent marginals). They illustrate the use of their
approach on small problems from portfolio management and queuing analysis, showing
that their uncertainty sets can provide constraint satisfaction guarantees equivalent to
classical models but are less conservative. To the best of our knowledge, their approach
have never been used in other studies or on other application cases.

1.3.0.2 Dimensionality reduction Ning and You (2018a,b) rely on Principal Com-
ponent Analysis (PCA) and Kernel Density Estimation (KDE) to build uncertainty sets
that directly integrate distributional information (correlation, asymmetry) available in
the historical data. They use PCA to identify correlation and KDE (as well as Robust
KDE) to define upper and lower bounds in each dimension of the principle components.
They illustrate the efficiency of this approach on different applications in robust model

26
predictive control, production batch scheduling and chemical process network design. In
the same framework, Zhang et al. (2022) use coupled Partial Least Square (PLS) or Kernel
PCA (KPCA) with Robust Kernel Density Estimation (RKDE) to construct uncertainty
sets.
Using dimensionality reduction and distribution estimation appears as a promising
approach to built uncertainty sets from historical data. However, PCA and similar tech-
niques must be used carefully as one has to set the number of principle components to
efficiently aggregate the information extracted from? the data. Moreover, PCA must be
applied on standardized data and ignore the difference in scale of uncertain parameters
which restricts its range of application to build uncertainty sets.

1.3.1 Gaussian Mixture Models

Campbell and How (2015) use the Dirichlet Process Mixture Model (DPMM) to define
uncertainty sets. This learning algorithm use historical data to build a Gaussian Mixture
Model that describe the distribution of uncertain parameters. The final uncertainty set
is the union of ellipsoidal uncertainty sets in which each component of the mixture model
is encapsulated in an ellipsoid. This approach leads to a SOCP robust counterpart and
have been applied by Zhao and You (2019) on a case study in energy production. Their
result shows that using the DPMM based uncertainty set succeed in reducing operating
costs against classical norm based uncertainty sets when distributions are multi-modal
and correlated and when the number of uncertain parameter is small (8). Ning and You
(2017) propose an alternative uncertainty set using DPMM. Their uncertainty set is the
union of several polyhedral (budget-based) uncertainty sets in which each component of
the mixture model is encapsulated in a bounded and convex polytope. They solve the
resulting robust optimization models using a cutting plane approach such as presented in
Section 1.1 in which several problems (one problem per component) are solved at each
iteration to identify violating uncertain realizations.

1.3.1.1 Support Vector Clustering Shang et al. (2017) introduced an uncertainty

set based on a one class Support Vector Clustering (SVC) algorithm. SVC is an unsuper-
vised version of the classical (One-Class) Support Vector Machines (SVM), a classification
machine learning algorithm that aims at defining the best boundary in order to partition a
space into two classes, based on a given set of (labeled) data points. Schölkopf et al. (2001)
propose to use SVC to describe the support of the distribution of (multi-dimensional)

27
random variables by defining cluster that contains a given portion of data. Shang et al.
(2017) build upon their approach and derive linear constraints to describe this cluster
as a bounded and convex polytope. When applied to the uncertain parameters of an
optimization problem, their method enables them to obtain new data-driven uncertainty
sets suitable for RO. The resulting uncertainty set efficiently captures correlations and
asymmetries in the distribution of uncertain parameters and the robust counterpart is
tractable (potentially large scale) when the initial problem is linear.
In the recent years, the SVC-based uncertainty set has attracted the interest of re-
searchers for various applications such as energy system optimization Shen et al. (2020),
resource allocation in a cellular network Wu et al. (2022), supply chain design Mohseni
and Pishvaee (2020), or model predictive control Shang and You (2019); Shang et al.
(2020). Several papers also extend the work of Shang et al. (2017) or propose alterna-
tive methods to build uncertainty sets based on support vectors. Mohseni and Pishvaee
(2020) extend the SVC-based uncertainty set with a Conservative Support Vector Cluster-
ing algorithm. Their approach provide higher satisfaction guarantees when only a small
data-set is available. In the same objective, Shang and You (2018) propose to reduce
the dimension of the uncertainty with PCA before using the SVC-based uncertainty set.
Han et al. (2021) and Lin et al. (2022) use Multiple Kernel Learning to improve the rep-
resentation accuracy of the support of uncertain parameters. The resulting uncertainty
sets are similar in their formulation to the classic SVC-based one but requires much more
variables and constraints, thus leading to very large robust counterparts and limiting the
range of application to problems with only a small number of uncertain parameters. Fi-
nally, Goerigk and Kurtz (2023) use One-Class Deep Support Vector Data Description to
formulate a SVC-based uncertainty set. This approach leads to non convex uncertainty
sets which can be used to model multimodal distributions more efficiently. However, the
learning process is extremely demanding and the resulting problem needs to be solved
with a decomposition approach (to tackle the non-convexity).

1.3.1.2 Sampling and Cutting Plane generation An other common approach in

DDRO is to consider the construction of uncertainty sets as a cutting plane generation
problem whose objective is to eliminate parts of the uncertain parameters space. In that
way, (Zhang et al. (2018b,a)) rely on copula theory to estimate the joint probability density
function (pdf) of uncertain parameters, then use Monte-Carlo sampling and define cutting
planes that eliminate a given portion of the samples. Zhang et al. (2018b) also compare

28
this approach with uncertainty sets based on the convex hull of the sampled realization.
All the uncertainty sets introduced in these papers are polyhedral uncertainty sets but
introduce a large number of cutting planes and can only be used to solve optimization
problems with a small (i.e. less than 10) number of uncertain parameters.
As an alternative to cutting plane generation, (Gumte et al. (2021b,a)) use Fuzzy
C-means clustering to define clusters in historical data and uniform sampling in these
clusters to define a discrete and non-convex uncertainty set.

2 Robust inventory management

This thesis focuses on the management of inventory in production systems subject to
uncertainties. In industrial systems, production decisions (e.g. quantity, scheduling) are
made in order to satisfy final customers demand and inventory decisions (e.g. placing
orders, quantity) are made to ensure the feasibility of the production plan.
It is well known that industrial systems are subject to various uncertainties (e.g. de-
mand, yield, processing times, etc) that affect their efficiency at each step of the produc-
tion process. In this context, robust optimization have been proposed as a method to
reduce the impact of uncertainties on production and inventory management decisions.
Due to the dynamic aspect of inventory management, most of the robust optimization
literature in this area have focused on multi-period problems. In this section we review
the main existing models and techniques used in the robust inventory management and
lot-sizing literature. To make this manuscript self supporting, we start with a brief intro-
duction to lot-sizing models in Section 2.1. We then present separately static and dynamic
approaches. The former, introduced in Section 2.2 corresponds to optimization problem
under uncertainty in which all the decisions are made prior to the uncertainty realization
(here and now decisions). Such problems can be formulated as the robust problem (1.1).
On the other hand, some of the decisions in dynamic formulations are adjustable, that is
they can be made once a portion of the uncertain parameters are realized. Section 2.3
presents applications of this principle to robust lot-sizing models.

2.1 Lot-sizing models

The Lot Sizing Problem (LSP) has been introduced by Wagner and Whitin (1958)
as a way to plan the production of goods (when and in which quantity) in order to

29
meet customers demand over a planning horizon in a deterministic context. Since then,
LSP models have been studied and extended to various types of production systems
(multi-echelon, multi-item, capacitated) and uncertainties (demand, lead-time, yield).
Since the last two decades, several papers proposed literature review on different aspects
of LSP models. Brahimi et al. (2006, 2017) provide a survey of the single-item LSP
model. Karimi and Fatemi (2003); Gicquel et al. (2008) propose reviews of capacitated
lot sizing models. Stadtler et al. (2015) propose a review of LSP models applied to
supply chain management and more recently, Metzker Soares et al. (2022) propose a
review of the different solution methods for LSP under uncertainty, including Stochastic
Programming, RO and Distributionally Robust Optimization. We refer interested readers
to these references for more information on LSP models.
The static and deterministic version of the single-item uncapacitated lot sizing prob-
lem aims at defining the setups (when to produce/order) and lot sizes (quantity pro-
duced/ordered) in a finite planning horizon of T periods in order to satisfy the customers
demands at minimum costs. The basic costs incurred consist in a setup cost s paid once
for each order, linear ordering costs v paid for each unit ordered and a linear holding cost
paid for each unit held in inventory at the end of a given period. We can formulate this
problem with the following mathematical model:

T
X
min (sxt + vqt + hIt ) (1.15)
t=1
t
X
s.t. It = qk − dk ∀t ∈ [T ] (1.16)
k=1

0 ≤ qt ≤ M xt ∀t ∈ [T ] (1.17)
xt ∈ {0, 1} ∀t ∈ [T ] (1.18)
It ≥ 0 ∀t ∈ [T ] (1.19)
(1.20)

in which for all t ∈ [T ], xt is a binary setup variable equal to 1 if production/ordering

occurs during period t, qt is the lot-size decision in period t, It is the inventory level at
the end of period t and M is a sufficiently large value (e.g. Tt=1 dt ).
P

In order to model the large variety of supply chains and production systems, this
basic model has been extended with various features in the literature. We presented
the single-echelon model in which a product is ordered or produced in order to satisfy

30
the customer demand with no intermediate sub-deliveries/assemblies. There also exists
a variety of multi-echelon models, in which several steps are necessary before meeting
the final customer demand. The same distinction is made for the number of item type
(single/multi-item).
Many robust inventory models in the literature include the possibility of a demand
shortage, where unmet demand can be either lost or delayed for a penalty cost. Specifi-
cally, lost sales corresponds to a case in which the unmet demand is penalized but can not
be satisfied in future periods. On the other hand when backordering or backlogging is al-
lowed, demand satisfaction may be delayed when the inventory is empty to future periods
one new units of product have been produced or received. Notice that some papers include
both type of penalties in a mixed approach to reflect practical applications (Nabil Absi
and Dauzère-Pérès (2011)). In some papers, both types of penalties are considered in a
mixed approach to reflect practical applications.
In practice, resources are often constrained by a capacity (e.g. maximum supplier
production, available production time, etc). When such constraints are applied, we talk
about capacitated models while uncapacitated models refer to problems with unconstrained
resources.
Finally, a wide variety of uncertainties can occur in supply chain and in production
systems. According to multiple review (e.g. Metzker Soares (2022)) the most studied
source of uncertainty is the demand. Other sources of uncertainty have also been studied
such as cost uncertainties (e.g. when production or ordering costs are uncertain) (Quezada
et al. (2019)), lead time uncertainties (Huang and Küçüyavuz (2008)), yield uncertainties
(Metzker Soares et al. (2022)) and production and setup times uncertainties (Taş et al.
(2019)).

2.2 Static robust inventory and production management

In this section, we provide an overview of existing works on robust inventory man-
agement static decision context, i.e. when decisions are made before the uncertainty is
realized. The first application of Robust Optimization in inventory management is Bert-
simas and Thiele (2004, 2006). The authors study a single-echelon, single-item LSP with
uncertain demand, as well as extensions to multi-echelon systems, with and without ca-
pacity on the production quantity or the inventory. They use a budget-based robust model
and derive a robust counterpart following the methodology introduced in Bertsimas and
Sim (2003, 2004). As a by-product, they propose a closed form policy in the case of a

31
single echelon problem without setup costs or with known ordering periods. Compared to
a dynamic programming (DP) solution, RO outperforms DP (better average and better
standard deviation) when the ratio backlog/holding is high enough and DP outperforms
RO in other cases. In the case of multi-echelon systems, RO outperform a myopic policy
(when DP is computationally intractable).
As a matter of example, the single-echelon single-item robust lot-sizing model of Bert-
simas and Thiele (2006) is given by:

T
X
min (sxt + vqt + yt ) (1.21)
t=1
t
( !)
X
s.t. yt ≥ max h qk − dk ∀t ∈ [T ] (1.22)
d∈Ut
k=1
t
( !)
X
yt ≥ max −b qk − dk ∀t ∈ [T ] (1.23)
d∈Ut
k=1

0 ≤ qt ≤ M x t ∀t ∈ [T ] (1.24)
xt ∈ {0, 1} ∀t ∈ [T ] (1.25)
yt ≥ 0 ∀t ∈ [T ] (1.26)

, where Ut = {z ∈ [0, 1]t |∥z∥1 ≤ Γt } and for all t ∈ [T ] the budgets Γt ∈ [0, t] are input
parameters that control the conservatism of the solution and where b is a backlogging
penalty.
Recently Hamed Mamani (2016) extends their work to more sophisticated uncertainty
sets (e.g. CLT based Bandi and Bertsimas (2012)) and Metzker Soares et al. (2022) used
the same methodology to LSP under yield uncertainty. These three papers derive closed
form policies that are useful for practitioners because they can be easily implemented on
large horizon without computational efforts. The use of such closed form policy is however
restricted to cases in which the ordering or production periods are known in advance as
they only compute the optimal ordering quantities.
A different solution approach was proposed in Bienstock and Ozbay (2008) for the
same problem without setup costs. They point out that considering constraints (1.22)
and constraints (1.23) can lead to sub-optimal robust decisions. In fact, this formulation
allows the inner maximization problems to fix the demand for a same period at two
different values in two distinct constraints. As a result, the worst holding/backlogging
cost yt can be overestimated and this can lead to over-conservative decisions. To overcome

32
this problem, Bienstock and Ozbay (2008) propose an iterative cutting plane algorithm in
which a single inner maximization problem is solved at each iteration. They illustrate the
efficiency of their approach by solving large instances of the problem (up to 300 periods).
Their works have been extended by Thorsen and Yao (2017) to the case with demand and
lead-time uncertainty, they use numerical experimentations to illustrate the conservatism
of the constraint dualization approach compared to the cutting plane approach.
More complex cases have also been studied in the context using a RO approach. Aouam
and Brahimi (2013) focus on an integrated production planning problem under demand
uncertainty in which the producer can decide to satisfy only a part of a client demand,
only incurring backlogging penalties on accepted orders. They also integrate the notion
of Work In Progress to model the congestion of production systems. Their experimental
results show that integrating order acceptance help to reduce congestion, increase the
fullfilment rate and increase companies profitability under uncertain demands. Varas
et al. (2014); Alem et al. (2018) propose a model to integrate lot-size and scheduling
decision in a robust context under demand uncertainty. Their model aims to define the
lot size of final product and schedule the usage of raw material with a time capacity for
the scheduling part. A nearly similar approach is presented in José Alem and Morabito
(2012) to combine lot sizing decisions and cutting stock problem. They face production
and holding costs, backlogging penalties but also integrate costs related to cutting stock
problem (e.g. trim loss). Finally, Agra et al. (2018) study a multi-product problem
directly inspired from an industrial case.
Most of the existing literature uses the classical budget based approach and does not
integrate more sophisticated structures of uncertainty sets such as presented in Section
1.2 and 1.3.

2.3 Dynamic robust inventory and production management

Dynamic solutions are applicable when a decision maker can adjust a part of its
decisions once the uncertainty is (at least partially) revealed. In inventory management,
such approaches allow to exploit the temporal dimensions of lot sizing problems and
provide less conservative solutions that classic RO formulations. One solution to correct
this lack of flexibility, Adjustable Robust Optimization (ARO) extends the RO approach to
a dynamic context by making some decisions adjustable to the realization of the uncertain
parameters. Multi-stage ARO suppose that the uncertainty is revealed over time and that
decisions must adjust to the current knowledge of the uncertainty. In the last two decades,

33
there has been many studies on ARO models, their computational complexity and how to
solve them. As this review focuses on papers centered on inventory or production systems
management that propose significant results on this topic, we refer interested readers to
the literature review presented in Daryalal et al. (2023) for a broader description of ARO.
In the multi-stage context implied by the structure of lot-sizing problems, ARO mod-
els often require to approximate the true adjustable problem to remain computationally
tractable. In that objective Ben-Tal et al. (2004) introduce the concept of affinely ad-
justable robust solutions, in which the adjustable decision variables are replaced by affine
functions of realizations of a subset of the uncertain parameteres. In the lot sizing prob-
lem for example, the quantity ordered/produced is expressed as a linear combination of
past demand realizations called linear decision rules. Let t ∈ [T ] and assume that we are
given the vector dt−1 = [d1 , . . . , dt−1 ] of realized demands up to period t − 1, the decision
variables qt , yt in period t then become:

t−1
qt0 qtk dk
X
qt (dt−1 ) = + (1.27)
k=1
t−1
yt (dt−1 ) = yt0 + ytk dk
X
(1.28)
k=1

Ben-Tal et al. (2004) apply this approach to solve a multi period, single item capacitated
LSP with multiple factories. They minimize the sum of all time-dependent production
costs in order to meet customers demand. They model the demand uncertainty with a
box uncertainty set and show the efficiency of their approach with experimental results.
The same methodology has been applied in Aharon et al. (2009) to multi-echelon sup-
ply chain problem with backordering and setup costs to study the impact of information
sharing. The same problem is studied in Melamed et al. (2016) under box and ellispoidal
uncertainty sets. See and Sim (2010) considered the same problem and propose to use
truncated linear decision rules. They provide models for the static case and both type of
decision rules under an ellispoidal uncertainty set but without any numerical results.
Finally, Solyalı et al. (2016) compare different formulations of static and dynamic
RLSP models. They compare both the quality of solutions and the tractability of models
under the formulations of Bertsimas and Thiele (2006), the adjustable policy of Ben-Tal
et al. (2004) and See and Sim (2010). They introduce a new formulation based on a
disaggregated model (i.e. facility location based LSP model). Bertsimas and Georghiou
(2015) compared the efficiency of affine, piecewise linear and piecewise constant decision

34
rules in a similar problem with supplier selection. At each period, an order is submitted
to a supplier with either instant delivery or a delay of period. Their numerical results
highlight that using several types of decision rules provides more flexibility in the nature
of adjustable decisions and outperforms other methods with a single type of decision rule.
Decision rules provide an efficient way to obtain a lower bound of the true ARO
problem by restricting the form of adjustable policies. In the opposite direction, Santos
et al. (2018) propose a lower bound for the ARO version of the lot sizing problem with
and without setup costs. They obtain a lower bound by solving a relaxed problem in
which the non-anticipativity constraints are relaxed and decisions are affine functions of
the demand on the whole planning horizon.
All studies considered up to now focuse on uncertain demand, which represents the vast
majority of the existing literature. However, recently other sources of uncertainty have
been studied in the context of ARO. Thevenin et al. (2022) propose an ARO approach to
the LSP with uncertain lead time (order cross-over), supplier selection and setup costs.
They use an exact solution method based on a constraint and column generation approach
(Zeng (2011)) and find through their numerical experiments that the problem is hard to
solve. To overcome this issue, they propose a fix and optimize heuristic that leads to a
significant reduction of the computation time.
Metzker Soares et al. (2023) study the dynamic lot sizing problem under yield uncer-
tainty and backorders via affinely adjustable robust policy. They show that the robust
counterpart based on affine decision rules is not tractable and compare different ways to
approximate it as a LP. More complex LSPs have also been studied via the ARO approach.
For example, Other works also include the relationship between multiple actors within a
supply chain. For instance, Ben-Tal et al. (2005) study a retailer-supplier commitment
contract problem in which the retailer makes prediction of its ordering quantities over the
planning horizon before any realization of the uncertain demands. These quantities may
be adjusted in a dynamic way after observing the demand, but induce additional penalties
when they deviate from its predictions (in addition to the LSP costs of setup, ordering,
holding and backorder). In their work, Ben-Tal et al. (2005) provide an ARO model for
this problem under classical box and ellipsoidal uncertainty sets. In this context, ARO
is an efficient way to obtain less conservative static decisions by reducing the amount of
items ordered during the planning horizon.

35
 Chapter 2

C ONTRIBUTIONS

In the previous chapter, we introduced the main robust optimization models used to
take decisions in supply chains and production systems. In this context, most papers focus
on uncertain demand while only two papers considers yield uncertainty (Metzker Soares
et al., 2022, 2023) and one paper focus on uncertain lead time and demand (Thorsen and
Yao, 2017). Moreover, in addition to the lack of diversity in the type of uncertainties,
most works have used classical box, budget-based or ellipsoidal uncertainty sets while only
three papers have used more complex uncertainty set (CLT based) and only one paper
has discussed the case of correlated uncertainty (Hamed Mamani, 2016). The use of more
sophisticated uncertainty sets integrating statistical information in robust models such
as described in Section 1.2 remains largely unexplored in the case of lot sizing models.
Moreover, one could exploit the amount of data collected by modern production systems
in order to propose a set of suitable uncertainty sets in order to handle the uncertainties
in a DDRO approach.
This gap in the literature may find its explanation in computational limitations. We
pointed out in Section 1.3 that despite their efficiency to model the behavior of uncertain
parameters and due to their more complex formulations, the use of data-driven uncertainty
sets comes at the expense of computational tractability. Such sets often lead to large scale
robust counterparts that are hard to solve in practice Zhang et al. (2018b); Han et al.
(2021); Bertsimas et al. (2017). This scalability issue is a limiting factor to the application
of DDRO in industrial problems often involving large scale mixed integer linear programs.
The rest of this chapter summarizes how the contributions detailed in Chapter 3, 4 and
5 attempt to address this issue.

Chapter 3: Approximate kernel learning uncertainty set In this first contribu-

tion we study the scalability of the data driven uncertainty set based on support vector
clustering. This uncertainty set has been applied successfully on several optimization
problems but suffers from scalability issues (see Chapter 1, paragraph 1.3.1.1). In fact,

37
the number of decision variables and constraints introduced in the RC is proportional to
both the number of uncertain parameters and the size of the data-set used to formulate
it. In practice, it results in large scale robust formulations that are practically intractable
when the initial problem is a MILP. In this chapter, we propose a two-phase method to
approximate the resulting uncertainty sets and overcome these tractability issues. The
main idea behind this approximation is to reduce the number of variables and constraints
introduced in the model by selecting a subset of historical data to construct an approxi-
mate uncertainty set. This subset is selected by an optimization process whose objective
is at first to maximize the similarity with the initial uncertainty set and then to minimize
the number of data point used in the formulation.
We evaluate our approximation method on three distinct well known optimization
problems, the robust binary knapsack problem, the robust lot-sizing problem as formu-
lated by Bertsimas and Thiele (2006) and a robust traveling salesman problem. Exper-
imental results show that the approximated uncertainty set leads to solutions that are
comparable to those obtained with the classic SVC-based uncertainty set on all problems
with a significant reduction of the computation time.

Chapter 4: A data driven robust optimization approach to supply planning.

In this chapter, we focus on a problem derived from a practical case encountered in
the aircraft industry. We consider a single assembly line that produces several distinct
item in order to satisfy the needs that are planned over a discrete and finite horizon.
The components are available in an external warehouse and the picking and delivery
of components to the assembly line are subcontracted to an external logistic provider
following the manufacturer orders. The problem consists in planning the delivery of
components on the assembly line and the quantity of each final item to produce at each
period of the planning horizon. The manufacturer follows the assemble to order rule,
which means that an item can not be produced earlier than its due date. In addition, the
time needed by the logistic provider to pick a set of components is uncertain and the total
picking time of the logistic provider during a given period is bounded by a capacity. If
the planning defined by the manufacturer exceeds this capacity in a given period, a part
of the components will not be delivered on time to the assembly line, possibly disrupting
the production plan and affecting its reliability. On the contrary, components delivered in
advance to the assembly line can reduce the efficiency of workers and the risk of incidents.
The objective is then to find a compromise between guaranteeing the feasibility of the

38
production plan and minimizing the stock held on the assembly line.
The manufacturer can track the state of each picking operation and the flow of compo-
nent from the warehouse to the assembly line leading to an historical data set of picking
and delivery time. We propose a robust mixed integer model based on a capacitated
lot sizing model to solve this problem. We compare the use of a classic budget based
uncertainty set and the SVC based uncertainty set to integrate picking time uncertainties
in the model. Building upon the contribution of Chapter 3, we use the approximation
of the SVC based uncertainty set to solve the data driven robust models in a reasonable
amount of time.
We conduct a numerical experiments on simulated data set inspired from the true
industrial case. The results show that the robust model can significantly reduce the
impact of picking time uncertainties on the assembly operations. In addition, the use of
the SVC-based uncertainty set leads to a lower level of component held on the border of
the assembly line while guaranteeing the reliability of the production plan.

Chapter 5: Data-driven Robust inventory management under uncertain de-

mand and lead time In this chapter we consider a more classic single item inventory
management problem under the two most common sources of uncertainty: demand and
lead time. From our previous studies, it appears that non-parametric models such as SVC
are particularly well suited to construct uncertainty sets when a large amount of data is
available. However, such large data-set are hard to obtain for demand or lead time uncer-
tainties as it would cover a very long time horizon during which the demand process can
evolve. In this context, we propose an alternative approach to build uncertainty sets from
a restricted set of historical data in which we assume that uncertain parameters follows a
specific (but unknown) stochastic process. We formulate an uncertainty set for demand
based on Autoregressive Moving Average (ARMA) process in order to provide a flexible
framework that can adapt to a wide range of demand processes. We conduced a set of
experiments for problems with lead time uncertainties, setup costs and autoregressive
demand. We compare the solutions obtained by a classic robust model using the budget
based uncertainty set, a sample average approximation model and our approach. Exper-
imental results show that our method results in less conservative solutions than a classic
budget based approach on average while maintaining a significant robustness advantage
over the sample average approach.
As mentioned in the previous chapter, using a dual reformulation approach to solve

39
this problem can be over-conservative as the same uncertain parameters appear in several
constraints simultaneously. Thus, we use a cutting plane approach such as described
in Thorsen and Yao (2017). In addition to a reduction of conservatism, this solution
approach provides several opportunities to improve uncertainty sets. In this context we
propose two improvements of the existing approach used to model lead time in robust
optimization: (i) We propose a discrete uncertainty set for lead times while the rest of the
literature have only considered continuous lead time models (Thorsen and Yao (2017);
Thevenin et al. (2022)), and (ii) we propose a simple approach to link uncertain lead
time and demand using linear regression in order to integrate potential correlation in the
variations of both types of uncertain parameters. In a second set of experiments, we
study the impact of uncertain lead time on the solution proposed by our three models.
Experimental results show that uncertain lead time reduces the gap of cost among all the
methods while preserving the conclusions drawn from the first set of experiments.

40
 Chapter 3

A PPROXIMATE K ERNEL L EARNING

U NCERTAINTY S ET

1 Introduction
As stated in Chapter 1, one of the main approaches to build data-driven uncertainty
sets is the use of support vectors. Since the introduction of this technique by Shang et al.
(2017) it has been applied in several applications (see Chapter 1). Unfortunately, all these
works agree on the lack of scalabilty of the method with respect to the number of uncertain
parameters. Indeed, so-called SVC-based uncertainty sets require the introduction of a
large number of new variables and constraints in the robust counterparts of optimization
models. These sets of variables and constraints grow proportionally to both the dimension
of the uncertainty and the number of data samples used to build them. In many practical
settings, in which we face a large optimization problem with a large data-set, this results
into intractable formulations. A first attempt to overcome this tractability issue was
proposed in Loger et al. (2022) to solve a robust inventory management problem under
SVC-based uncertainty set. As an example, Figure 3.1b shows the increase in CPLEX
CPU time to solve the robust 0-1 Knapsack problem with 30 items under uncertain weights
with respect to the number of data samples.
This chapter proposes a method to reduce the number of additional variables and
constraints associated with an SVC-based uncertainty set. More precisely, we develop a
method to approximate the SVC-based uncertainty set introduced in Shang et al. (2017) in
order to overcome the scalability issue that was pointed out. The solutions computed with
the resulting approximate uncertainty set are identical or comparable to those obtained
with the SVC-based uncertainty set with a significant reduction of the computational
burden.
The rest of this chapter is organized as follows: Section 2 describes the uncertainty

41
 (a)

(b)

Figure 3.1 – Performance of the SVC-based uncertainty set on a robust 0-1 Knapsack
problem with uncertain weights and 30 items. Top: Trade-off between solution quality
and feasibility using the SVC-based (USVC ) and the budget-based (UΓ ) uncertainty sets.
Bottom: Evolution of CPLEX CPU time in seconds using the SVC-based uncertainty set
for different sizes of data sets.

42
set of Shang et al. (2017) and provides some explanations on its scalability issues. We
introduce our approximation procedure in Section 3 and present experimental results in
Section 4. The latter first focuses on the ability of the approximation to define a set
that resembles the original SVC-based uncertainty set, then provides evidence of the
optimization and computational performance on three different optimization problems.

2 Support vector clustering-based uncertainty set

SVC is a well-known Machine Learning algorithm, introduced by Shang et al. (2017),
that describes the support of a multi-dimensional random variable with a cluster of min-
imal volume containing a given portion of historical realizations. Among the various ad-
vantages displayed by SVC, its usage as a ML tool to define uncertainty sets is motivated
by the fact that it naturally incorporates information on mixed moments of distributions
such as the covariance of random variables. In addition, depending on the chosen Kernel
function, it does not necessitate any complex tuning of hyper-parameters. To make this
chapter self-contained, we briefly describe below the approach developed in Shang et al.
(2017).
n o
Assume that we dispose of a set D = u(1) , . . . , u(N ) of N data samples in an m-
dimensional space. The principle of SVC is to use a non-linear mapping ϕ(u) : Rm → Rk
to a high-dimensional feature space F and to seek the smallest sphere that encloses at
least (1−ν)N data samples, where the parameter ν ∈]0, 1[ is a lower bound on the portion
of data samples considered as outliers. The general idea is to express this optimization
problem as a Lagrangian function that is then minimized using the KKT conditions.
Shang et al. (2017) apply the so-called kernel trick with a suitable piecewise linear function
called WGIK (Weighted General Intersection Kernel) to obtain a quadratic program (QP)
that computes the best sphere with respect to D. The different step used to obtain the
QP and its interpretation are given in Appendix A.
The QP solution defines a subset S of data points, called support vectors (SV), that
lie outside of the boundary of the sphere. A subset B ⊆ S of these points are exactly on
the boundary of the sphere and are referred to as boundary support vectors (BSV). To
summarize, those two sets are such that B ⊆ S ⊆ D and |S \ B| ≤ N ν ≤ |S|. In the
robust optimization framework, the parameter ν is diverted from its original purpose to
control the size of the uncertainty set and consequently the level of conservatism of the
model. That is, selecting a larger value of ν increases the number of SV and decreases

43
the conservatism. Figure 3.2 represents these sets for 1000 samples following a bivariate
gaussian distribution when ν is set to 0.15.

Figure 3.2 – Representation of S and B for a bivariate gaussian distribution when ν = 0.15

From now on, we slightly abuse the notation and sometimes refer to the data point
u using its index i, e.g., writing i ∈ S instead of u(i) ∈ S, when the meaning is clear
(i)

from the context. For each point i ∈ D, the QP also computes a value αi ∈ [0, 1/νN ],
where αi = 0 iff i ∈ D \ S, αi = 1/νN iff i ∈ S \ B and 0 < αi < 1/νN if i ∈ B. Let
α = [α1 , . . . , αN ] be the vector of obtained coefficients, Shang et al. (2017) use α to define
the following data-driven uncertainty set:
( )
(i)
X
USVC = u αi ∥Q(u − u )∥1 ≤ θ (3.1)
i∈S

P ′

where ∥·∥1 stands for the ℓ1 -norm, Q is a weighting matrix and θ = min
′ i∈S αi ∥Q(u(i ) − u(i) )∥1 .
i ∈B
|S|
For simplicity, we also define for all ᾱ ∈ R+ and u ∈ Rm
+ the function:

ᾱi ∥Q(u − u(i) )∥1

X
θ(ᾱ, u) = (3.2)
i∈S

which measures the weighted distance of u to all points in S used in definition (3.1) of
USVC . Note that USVC = {u |θ(α, u) ≤ θ }. In particular, there exists a data point i′ ∈ B
′
such that θ(α, u(i ) ) = θ, where θ is the threshold used to define the boundaries of USVC .

44
In what follows, we extend definition (3.1) and denote
( )
(i)
X
U ᾱ = u ᾱi ∥Q(u − u )∥1 ≤ θ = {u |θ(ᾱ, u) ≤ θ } (3.3)
i∈S

for any vector ᾱ ∈ Rm

+ . In particular when ᾱ = α, we have USVC = Uα. An appealing
property of this class of uncertainty sets is that one can easily reformulate U ᾱ as a
polyhedron using auxiliary variables V = [v1 , . . . , v|S| ]:
 




∃vi ∈ Rm + , i ∈ S s.t




U ᾱ = u T
i∈S ᾱi vi 1 ≤ θ (3.4)
P
 

−vi ≤ Q(u − u(i) ) ≤ vi , i ∈ S

 


which is bounded and nonempty for 0 < ν < 1 (Shang et al. (2017)). Based on formulation
(3.4), the left-hand side of the robust linear constraint

max aT x ≤ b (3.5)
a∈USVC

with a ∈ USVC ⊆ Rm and b ∈ R, is thus equivalent to the following LP:

max aT x
αi 1T vi ≤ θ
X
s.t.
i∈S

−vi ≤ Q(a − u(i) ) ≤ vi ∀i ∈ S

which is feasible and bounded since USVC is nonempty and bounded for 0 < ν < 1. By
strong duality, the dual of this problem

(µi − λi )T Qu(i) + ηθ
X
min (3.6)
i∈S
X
s.t. Q(λi − µi ) + x = 0 (3.7)
i∈S

λi + µi = ηαi 1 ∀i ∈ S (3.8)
λi , µi ∈ Rm
+ ∀i ∈ S (3.9)
η ≥0 (3.10)

is also feasible and bounded and their optimal value coincide. Therefore, the robust

45
counterpart of constraint (3.5) is given by

(µi − λi )T Qu(i) + ηθ ≤ b
X
(3.11)
i∈S

with the additional constraints (3.7)-(3.10).

As the formulation (3.11), (3.7)-(3.10) suggests, the size of the robust counterpart
depends upon the number of support vectors. Specifically, consider a robust constraint
j that contains mj uncertain parameters taking their values in an SVC uncertainty set
defined with support vectors Sj . The robust formulation then adds 2|Sj | · mj + 1 new
variables and (|Sj | + 1) · mj new constraints to the original one. Since the number |Sj |
of support vectors is closely related to νN , the size of the robust counterpart is greatly
influenced by both the regularization parameter ν and the size of the original dataset D.
In practice, increasing one of these values may greatly affect the tractability of the original
program and become a significant limitation to the application of this framework. This
phenomenon is particularly problematic in the context of big data where one would need
to take advantage of information contained in a high volume of collected data. Moreover,
it is well known that the SVC algorithm suffers from the curse of dimensionality, as
stated by Keogh and Mueen (2017): « For machine learning problems, a small increase
in dimensionality generally requires a large increase in the numerosity of the data, in
order to keep the same level of performance for regression, clustering, etc. ». In our case,
increasing the number of uncertain parameters mj requires more data samples for the
SVC algorithm to be efficient, leading in turn to a large robust counterpart unlikely to be
solvable in a reasonable time.

3 Approximation of the SVC-based uncertainty set

To overcome the tractability issue pointed out in the previous section, we propose a
two-step method to reduce the size of the robust program derived from the SVC-based
uncertainty set. The goal is to curtail the number of support vectors that define USVC ,
while keeping a good representation of the original cluster obtained via the SVC.
Assume that, for a given ν, coefficients α and the corresponding SVC-based uncer-
tainty set USVC are computed following the methodology of Shang et al. (2017). In a nut-
shell, the procedure consists of computing a vector of modified weights α̂ = [α̂1 , . . . , α̂N ]
that is used to derive an approximation of the original uncertainty set. The algorithm

46
takes an integer K as an input parameter, which sets the maximum number of strictly
positive coordinates of α̂, thereby limiting the number of SV that shape the (modified)
uncertainty set. This in turn reduces the size of the resulting LP (3.6)-(3.10) and its
robust counterpart, since we only need to include dual variables for SV u(i) with a strictly
positive coefficient α̂i in their formulation. After selecting an initial subset of SV, the al-
gorithm runs in two phases: (1) optimize the weights α̂ of the selected points in order to
minimize a given distance between U α̂ and USVC , then (2) minimize the number of points
selected without deteriorating the distance obtained during the first step. We detail the
procedure below.

3.1 Initial support vectors selection by K-medoid clustering

As a strategy to reduce the size of the set S, the algorithm starts by choosing a first
subset Sb0 of support vectors that are a priori good representatives of S. By default, we
include all the boundary support vectors B in Sb0 , as well as K additional support vectors
selected from S \ B. The general approach consists of partitioning S \ B into K distinct
subsets of data points sharing similar features and aggregate the information contained
in each subset into a single data point. The underlying idea is that close samples that lie
strictly outside USVC have similar contributions to cluster boundaries.
We perform the partition using the so-called K-medoid clustering algorithm described
in (Schubert and Rousseeuw, 2019), also known as partitioning around medoids (PAM,
see Kaufman and Rousseeuw (1990)). This procedure assigns each sample in S to exactly
one cluster in order to minimize its distance to the center of this group. One appealing
property of the K-medoid clustering compared to the more well-known K-means cluster-
ing is that it forces the center of each cluster to be an actual data point. In our case, this
ensures that any selected point belongs to S \ B and therefore lies outside of the initial
uncertainty set USVC . We compute the distances between two points u, v ∈ Rm using the
function
d(u, v) = ∥u − v∥2 (3.12)

where ∥·∥2 represents the ℓ2 -norm. At the end of the K-medoids procedure, we obtain K
clusters along with their respective centers. The union of B with these centers form the
initial set Sb0 , as represented in Figure 3.3b. In our numerical experiments, we rely on an
existing implementation of the algorithm available in julia 1 .
1. https://juliastats.org/Clustering.jl/dev/index.html

47
3.2 Procedure OPT -α
Recall from the definition (3.1) that we use the quantity θ(α, u(i) ) = i∈S αi ∥Q(u −
P

u(i) )∥1 to determine whether a point u ∈ Rm + belongs to USVC . Following the reduction of
the original set of support vectors into the smaller subset Sb0 , our goal is now to define an
alternative uncertainty set U α̂ where α̂i > 0 iff i ∈ Sb0 . A simple modification would be to
set α̂i = αi 1{i∈Sb0 } , but this is likely to distort the boundaries of the resulting uncertainty
set. Instead, we propose to compute modified weights α̂ such that the uncertainty set U α̂
approximate USVC in a satisfactory fashion.
Our goal is then to limit as much as possible the size of the subset (USVC \ U α̂) (U α̂ \ USVC ).
S

Since a direct calculation of the latter would clearly be too cumbersome, we instead use
the initial sample of data points to measure the quality of the approximation U α̂.
Specifically, we rely on the empirical cumulated absolute deviation

|θ(α̂, u(i) ) − θ(α, u(i) )|

X
∆= (3.13)
i∈D

Notice that all points i ∈ S that are discarded in Sb0 are not candidate to serve as
support vectors. Hence, we artificially set their modified weights α̂i to 0. We use ∆ as
a mean to evaluate the quality of the obtained approximation. We seek values α̂ that
minimize this empirical distance, which amounts to solving the optimization problem
 
X X 
(i′ ) (i′ )
− u(i) )∥1 − − u(i) )∥1 
X
min  αi ∥Q(u α̂i ∥Q(u (3.14)
α̂
i′ ∈D i∈S i∈S

It is straightforward that (3.14) is equivalent to the following linear program:

X
min ∆i′ (3.15)
i′ ∈D
′
(αi − α̂i )∥Q(u(i ) − u(i) )∥1 ≤ ∆i′ ∀i′ ∈ D
X
s.t. (3.16)
i∈S
′
(α̂i − αi )∥Q(u(i ) − u(i) )∥1 ≤ ∆i′ ∀i′ ∈ D
X
(OPT -α) (3.17)
i∈S

α̂i = 0 ∀i ∈ S \ Sb0 (3.18)

α̂i ≥ 0 ∀i ∈ Sb0 (3.19)
∆i′ ≥ 0 ∀i′ ∈ D (3.20)

48
where constraints (3.16) and (3.17) linearize the absolute distance used in (3.14) and the
variables ∆i′ are the absolute deviations described in (3.13) evaluated for a given point
i′ ∈ D. The solution to this LP provides the optimal weights α̂0 that we can then use
to define the modified uncertainty set UKM = U α̂0 . This first approximation is illustrated
on Figure 3.3b where the orange triangles represent the medoids that constitute Sb0 .

3.3 Reduction of the number of support vectors

The objective of the second phase is to further reduce the number of support vectors
retained in Sb0 while keeping a good quality of the first approximation. At this stage, we
relax the constraint that a support vector must belong to the initial set of medoids or to
B. Instead, we consider any subset of points Sb ⊂ S as potential support vector for the
approximation. Let ∆∗ be the value of the optimal solution to (OPT -α) computed in the
first phase. We impose that the solution derived from this second phase: (1) selects at
b ≤ |Sb |, and (2) defines new weights α̂∗ such that
most K + |B| support vectors, i.e. |S| 0

α̂i > 0 iff i ∈ S and the resulting empirical distance between USVC and U α̂∗ is at most
∗ b

∆∗ . We introduce the set of binary variables


 1 if u(i) ∈ Sb
yi =
 0 otherwise

and formulate the above problem with the following MIP:

X
min yi (3.21)
i∈S
′
(αi − α̂i )∥Q(u(i ) − u(i) )∥1 ≤ ∆i′ ∀i′ ∈ D
X
s.t. (3.22)
i∈S
′
(α̂i − αi )∥Q(u(i ) − u(i) )∥1 ≤ ∆i′ ∀i′ ∈ D
X
(3.23)
i∈S

∆i′ ≤ ∆∗0
X
(3.24)
i′ ∈D

α̂i ≤ yi ∀i ∈ S (3.25)
α̂i ≥ 0 ∀i ∈ S (3.26)
∆i′ ≥ 0 ∀i′ ∈ D (3.27)
yi ∈ {0, 1} ∀i ∈ S (3.28)

49
where (3.21) minimizes the number of SV retained in S. b Constraints (3.22)-(3.24) ensure

that the empirical distance achieved by the selected support vectors with their modified
weights is no greater than ∆∗ , i.e,. that the quality of the first approximation is preserved.
Finally, constraints (3.25) ensure that the weight α̂i associated to a sample point u(i) is
positive only if the variable yi is equal to one, that is if it belongs to S.
b

The solution to the MIP above defines a new subset of support vectors Sb that we
can use to approximate USVC . Notice that the optimal weights α̂0 to (OPT -α) directly
define a feasible solution of value K + |B| by letting yi = 1 for all i ∈ Sb0 ∪ B and yi = 0
otherwise. It follows that the approximation UKM can be used as an initial solution to the
MIP (3.21)-(3.28). The corresponding initial value is thus an upper bound, which ensures
that the final set of Sb of SV satisfies |S|
b ≤ |Sb |. Notice that at this point, the value of
0

the modified weights α̂ may not be optimal since this MIP primarily focuses on reducing
the size of S.
b However, running the program (OPT -α) on this new subset Sb allows us

to compute optimized weights α̂∗ in order to minimize the distance ∆ defined by (3.13).
Let UAC = U α̂∗ be the resulting uncertainty set. Figure 3.3c represents the uncertainty
set UAC and the set Sb ⊆ S.

3.4 Summary of the approach

We summarize the procedure to construct this new uncertainty set U α̂ in three steps:

Initialization Apply the method of Shang et al. (2017) to obtain the weights α and the
associated sets S and B.

Step 1 Apply the K-medoid algorithm to obtain an initial subset Sb0 ⊆ S, then solve
(OPT -α) to derive the lower bound ∆∗ and collect the modified the weights α̂0 to
define UKM .

Step 2 Solve the MIP (3.21) - (3.28) to obtain the final subset S,
b then use (OPT -α) to

refine the coefficients α̂∗ and obtain UAC .

One can then follow the method described in Section 2 to obtain the robust counterpart
of a constraint with the new uncertainty set UAC = U α̂∗ . Figure 3.3 illustrates the
evolution of the set of SV and the boundaries of the resulting uncertainty set throughout
the different step of the procedure on a simple two-dimensional example, in which the two
uncertain parameters follow a bivariate normal distribution.

50
 (a) Initial cluster

(b) First approximation

(c) Final approximation

Figure 3.3 – Steps of uncertainty set construction

51
4 Experimental results
In this section we present several numerical experiments to evaluate the performance
of the approximation methods introduced in this chapter.
The first part compares the approximation UAC derived from the procedure intro-
duced in the previous section with the original SVC-based uncertainty set USVC . The
goal is to produce numerical evidence that the proposed optimization procedures lead to
a satisfactory approximation in the sense that USVC and UAC define comparable subsets
of Rm . In the second part, we apply the resulting uncertainty set on three distinct robust
optimization problems. We find consequent performance improvements without signifi-
cantly deteriorating the solution quality. Specifically, we first consider a robust knapsack
problem with uncertain weights similar to the one studied in Bertsimas et al. (2004), a
robust Traveling Salesman Problem (TSP) with travel time uncertainties, and the robust
inventory management problem of Chapter 3 in which the demand is an autocorrelated
stochastic process. For each problem we examine (1) the quality of the solutions obtained
by the original SVC-based uncertainty set USVC , and (2) the reduction of the computa-
tional burden incurred to solve the problem to optimality, compared to USVC . As the MIP
(3.21) - (3.28) is hard to solve to optimality, in what follows we restrict the time allowed
to solve it to 1 hour. All models are implemented in Julia 1.9.0 with JuMP and solved
with CPLEX 20.1 using a single thread on a Linux, Ubuntu 20.04.2 LTS, equipped with
an Intel Xeon Gold 6230 clocked at 2.10GHz.

4.1 Uncertainty set comparison

In this section, we compare the portion of the distribution of random parameters
covered by USVC and the approximate uncertainty sets obtained after the second phase of
the procedure. Our test bench covers different numbers of uncertain parameters and uses
two types of distributions (Gamma and Gaussian). For each experiment, we generate a
set D of N data points that serves as input data to define the considered uncertainty sets
(USVC and UAC ) and a distinct test set T of 10000 data points to evaluate the out-of-
sample error. Finally, we evaluate the portion of T contained in USVC and UAC . To avoid
special cases and strengthen the conclusions, we repeat the above procedure 10 times.
The cumulative out-of-sample error (Type I and II) is given as a percentage of |T |. In
what follows, the Type I error refers to the percentage of data points in T that belongs
to UAC but not to USVC and the Type II error refers to the percentage of data points in

52
Figure 3.4 – Evolution of the out-of-sample cumulative error as a function of the ratio
K/m for a multivariate Gamma distribution when N = 1000

T that belongs to USVC but not to UAC .

Figure 3.4 represents the cumulative error computed for UAC on the test sample T for
various values of m (number of uncertain parameters) and several ratios K/m. We repre-
sent the average error for various values of ν ∈ {0.15, 0.25, 0.35, 0.45, 0.55, 0.60, 0.70, 0.75,
0.80, 0.85, 0.90, 0.95, 0.99}. It appears that the parameter K offers an efficient control over
the quality of the approximation.
Figure 3.5 displays the average percentage of total error obtained as a function of ν.
We find that this percentage increases as ν increases from 0.15 to 0.8 but decreases on
[0.8, 0.99].
Remarks:

1. As mentioned above, the MIP (3.21)-(3.28) proved to be difficult to handle by

commercial solvers like CPLEX or Gurobi. In our experiments we let 1 hour to the
solver to improve upon the initial solution. In comparison to this time limit, the
computation time required by Step 1 (i.e. the K-medoid algorithm and solving LP
(OPT − α)) is negligible. On all the approximate uncertainty sets computed for
this set of experiments, it took on average less than 1 second and always less than
3 seconds to complete Step 1.

2. We conducted other experiments with larger data-sets (up to N = 5000) and there
was no significant difference in the quality of the approximation UAC .

53
Figure 3.5 – Evolution of the cumulative error as a function of ν for a multivariate Gamma
distribution when N = 1000, K/m = 3

4.2 Application to optimization problems

The robust problems that we consider in this section are expressed as MIP of the form

max c⊺ x (3.29)
s.t. u⊺j x ≤ bj ∀uj ∈ Uj (D) (3.30)
x ∈ {0, 1}m (3.31)

where uj is the vector of uncertain parameters of the jth constraint, which are re-
stricted to belong to the uncertainty set Uj (D). In our computational experiments, we
distinguish four possible implementations of Uj (D), all derived from the data available
D = {u(1) , . . . , u(N ) }. Each type of uncertainty set results in a different robust counter-
part for each of the problems under study. In the following analysis, we denote by zSVC
and zAC the final objective function value attained with the formulation based on the
uncertainty set USVC and UAC , respectively. Similarly, we let τSVC and τAC be the needed
absolute CPLEX CPU-time by both formulations to converge. We compare the respec-
tive performance of the two uncertainty sets based on the relative objective value function
ρ = (zSVC −zAC )/zSVC , while the relative CPLEX CPU-time is given by (τAC /τSVC )·100%.
In the lot sizing problem and in the TSP, the uncertainty affects the objective function,
hence zSVC and zAC are evaluated a posteriori on an independent test set of 10000 new
data samples. Note that for the 0-1 Knapsack problem, some realizations of the uncertain
parameters may lead to the constraint violation. In that case, we also use a distinct test
set of 10000 new samples to compute the achieved fraction of feasible solutions φSVC and

54
φAC using USVC and UAC , respectively. We then express the relative constraint satisfaction
as the ratio (φSVC − φAC )/φSVC .

4.2.1 Robust 0-1 Knapsack Problem

The first problem considered in this section is the robust knapsack problem with
uncertain weights.
This problem has been introduced in Bertsimas et al. (2004) and can be expressed
with the following IP:

max c⊺ x (3.32)
s.t. w⊺ x ≤ W ∀w ∈ U(D) (3.33)
x ∈ {0, 1}m (3.34)

where W is the total capacity and w ∈ Rm is the vector of weights of the items considered,
which is restricted to the uncertainty set U(D).
We generate instances of different size (i.e., number of items). Random weights
are generated using multivariate gamma distributions with random correlations drawn
uniformly in [−1, 1], shape and scale parameters drawn uniformly on [10,13]. As in
Bertsimas et al. (2004), costs are randomly generated in [17, 77]. We considered in-
stances of size m ∈ {10, 20, 30, 40} and generate 10 instances of each size. Note that
we solve 8 distinct models on each instance, each of them corresponding to a specific
value. ν ∈ {0.15, 0.25, 0.35, 0.55, 0.75, 0.85, 0.95, 0.99}.Figure 3.6 represents, respectively,
the difference of objective value ρ and constraint satisfaction between solutions to models
based on USVC and UAC for different number of uncertain parameters and different value
of the ratio K/m. Regarding the differences in terms of constraint satisfaction, we ob-
serve that our approximation leads to a higher level of constraint satisfaction than the
SVC-based model although the difference becomes less and less significant as the ratio
K/m increases. On the other hand, our approximation displays smaller values of objective
function in most cases, where again the difference decreases with the ratio K/m.
Figure 3.7 represents the absolute and relative solving time of UAC for different values
of m and of the K/m ratio. We observe that using this approximate uncertainty set leads
to a reduction of the solving time by more than 50% on all models. This phenomenon
diminishes as the number of uncertain parameters increases: Since the number of SV used
to define USVC only depends on ν and the number of data samples while our approximation

55
Figure 3.6 – Relative constraint satisfaction (left) and objective function (right) difference
in percent for the knapsack problem

Figure 3.7 – Absolute (left) and relative (right) solving time of UAC as a function of m for
the robust knapsack problem with different values of the ratio K/m

56
Figure 3.8 – Trade-off between objective value and constraint satisfaction of USVC and
UAC for different values of K/m

scales with the dimension m, the relative advantage of the latter is less significant when
these two values compare less favorably.
Figure 3.8 represents the constraint satisfaction level in percent as a function of
the objective value of the solution obtained with USVC and UAC . The considered in-
stance has 30 items and the different solutions correspond to different values of ν ∈
{0.15, 0.25, 0.35, 0.45, 0.55, 0.75, 0.85, 0.95}. This example illustrates the typical behavior
observed for solutions obtained using USVC and UAC in our experiments. The majority of
the marked points (corresponding to a specific value of ν) display comparable performance
between the original robust solution and the approximated one, even when K/m = 1. In
particular for K/m = 3, both the objective value and the probability of constraint satis-
faction are equivalent when using UAC instead of USVC .
Figure 3.9 represent the number of support vectors |S|b used in the approximation

UAC as a function of m and different values of the K/m ratio. When the K/m ratio is
less than or equal to 2, we observe that the number of support vectors is less than m
for most instances. When K/m = 3, the average number of support vectors is less than
1.5m. Finally, when the ratio K/m = 4, the number of support vectors is less than 4m
but the second phase does not seem to reduce that value significantly within the allowed
computation time of 1 hour.

4.2.2 Robust Lot-sizing Problem

In this section we consider an uncapacitated lot-sizing problem (ULSP) similar to the

one addressed in Bertsimas and Thiele (2006). The goal is to define the ordering quantities

57
Figure 3.9 – Number of support vectors |S|
b used in U
AC as a function of m and different
values of K/m

of a single item over a finite and discrete planning horizon of T periods. In each period
t = 1, . . . , T , we assume an uncertain demand dt that can be satisfied immediately using
units held in the inventory or (fully or partially) backlogged to be served in a subsequent
period. Four distinct types of costs are incurred: (1) a per-unit ordering cost c, (2) a
fixed ordering cost s that is paid in every period when an order is placed (regardless of
the quantity), (3) a per-unit, per-period holding cost h that applies to units physically
held in the inventory to satisfy future demands, and (4) a per-unit, per-period backlogging
penalty cost b for every unmet unit of demand. The objective is to minimize the total
cost incurred by the system of the planning horizon [T ].

Bertsimas and Thiele (2006) model this problem with the following MIP:

T
X
min (cqt + Sxt + yt ) (3.35)
t=1
t
!
⊺
X
s.t. y t ≥ h x0 + ui − min dt 1 ∀t ∈ [T ] (3.36)
d∈Ut (D)
i=1
t
!
⊺
X
yt ≥ −b x0 + ui − max d1 ∀t ∈ [T ] (3.37)
dt ∈Ut (D)
i=1

0 ≤ qt ≤ M x t ∀t ∈ [T ] (3.38)
xt ∈ {0, 1} ∀t ∈ [T ] (3.39)
yt ≥ 0 ∀t ∈ [T ] (3.40)

58
where qt is the quantity of items ordered in period t, xt is a binary variable equal to
1 if qt > 0, 0 otherwise, yt correspond to the holding/backlogging cost of period t ,
and M is a large constant that can be set to Tt=1 dt w.l.o.g.. In constraints (3.36)-(3.37),
P

Ut (D) denotes one of the data-driven uncertainty set that describes the uncertain demand
parameters dt = [d1 , . . . , dt ] that contains the demand values over periods 1 to t.
In our experiments, we consider a horizon length of T = 14 and let S = 3000, c = 1
and h = 4. We consider random demand scenarios inspired from the ones introduced in
Ben-Tal et al. (2004) that we generate as follows: for all period t ∈ [T ], we first define a
seasonal nominal demand
!!
1 2π(t − 1)
d¯t = 1000 · 1 + sin
2 7

and use it to generate an autocorrelated final demand process (dt )t∈[T ] with autocorrelation
factors a1 , . . . , ak and error term dˆt :

k
dt = d¯t + aj dt−j + dˆt .
X

j=1

The perturbations dˆt are independent and identically distributed according to either
a Gamma or a Normal distribution and satisfy for all t ∈ [T ] the property P (−0.2 ≤
dˆt ≤ 0.2) ≥ 0.9. For each distribution, we generate 5 different data sets and 3 different
instances with varying ratio σ = hb ∈ {0.5, 1, 3}. We thus obtain 30 instances of the
ULSP with stochastic demand. As in the previous experiments, each instance is solved
for different values of ν ∈ {0.15, 0.25, 0.35, 0.55, 0.75, 0.85, 0.95, 0.99} and each solution is
evaluated on a data set containing 10, 000 data points based on the ratio ρ.
Figures 3.10 show the average and standard deviation of the difference in objective
value for several values of the ratios K/m and σ. Notice that we present aggregated
results for the two considered types of distribution and all possible values of ν since
our experiments did not reveal any significant differences when these parameters vary.
For all values of σ, both the average and standard deviation (across all instances) of ρ
decrease as K/m increases. Most of the solutions computed based on UAC are worse than
those obtained with USVC , although the average gap and standard deviation between the
two diminishes quickly with K/m. Again, we observe that a ratio K/m = 3 provides
a very good approximation. We also observe that the difference in objective value is
less consistent than the one obtained for the Robust 0-1 Knapsack problem. This is in

59
Figure 3.10 – Average (left) and standard deviation (right) of the relative difference be-
tween the objective values zSVC and zAC for the robust lot sizing problem.

Figure 3.11 – Solving time of UAC in percent of the solving time of USVC as a function of
K/m for different values of σ (Robust lot sizing problem)

contrast with a problem using binary or integer variables, for which a slight change of
the boundaries of the uncertainty set is less likely to affect the feasibility of the optimal
solution.
Figure 3.11 shows the solving time of UAC in percent of the solving time of USVC for
different values of K/m and different ratios σ. We again observe that the approximate
uncertainty set succeeds in reducing the time needed to solve the robust counterpart as
this ratio never exceeds 4%. For the SVC-based model, the median of CPLEX CPU-time
is above 1000 seconds and goes up to 10, 000 seconds in the slowest cases. In comparison, it
is always shorter than 10 seconds when the robust counterpart uses UAC even for K/m = 4.
In addition, we observe a slightly higher computation time for both uncertainty sets when
the ratio σ = hb is equal to 0.5.

60
4.2.3 Robust Traveling Salesman Problem

In this section, we consider a robust traveling salesman problem with uncertain travel
costs. Given a set of m cities, the objective is to compute a tour of minimal cost that
visits each city once and returns to the starting city. We consider the case where some of
the costs cij incurred when traveling from city i to city j are uncertain. We denote by A
and à the sets of arcs with known and uncertain traveling costs, respectively.
For all i, j ∈ [m], let xij be a binary variable equal to 1 if the arc (i, j) is used in
the solution and ui ∈ [m] be an integer variable that gives the position of city i in the
corresponding tour. We rely on the well-known Miller-Tucker-Zemlin (MTZ) formulation
to express the problem as a MIP:

min Z (3.41)
X X
s.t. xij cij + max xij cij ≤ Z (3.42)
c∈U (D)
(i,j)∈A (i,j)∈Ã
X
xij = 1 ∀j ∈ [m] (3.43)
i∈[m]
X
xij = 1 ∀i ∈ [m] (3.44)
j∈[m]

ui − uj + mxij + 1 ≤ m 2 ≤ i ̸= j ≤ m (3.45)
1 ≤ ui ≤ m ∀i ∈ [m] (3.46)
xij ∈ {0, 1} ∀i, j ∈ [m] (3.47)

where U(D) is one of the data-driven uncertainty sets built from historical traveling costs.
We generate instances with m ∈ {15, 20, 25} cities and different partitions A and
Ã. The coordinates of cities are generated randomly and drawn uniformly in the square
[1, 50]2 . The traveling time on edges belonging to A and the average traveling cost for
edges belonging to à are set to the Euclidean distance between the two extremities. To
select the set of edges that belong to Ã, we used the following procedure:
1. Select n cities randomly in [m]
2. For each of the n selected cities compute the set C of the 3 closest cities in terms of
Euclidean distance.
3. Add edges (j, i) and (i, j) to Ã.
In this set of experiments, we used values of n ∈ {5, 7}. Traveling costs are such that
cij = c̄ij +δij ∀(i, j) ∈ Ã where c̄ij is a fixed traveling cost and δij is a random perturbation

61
following a Gamma (for half of the instances) or a Normal distribution (for the other half).
For each value of m, n and each type of distribution, we generate 2 instances for a total of
24 different instances. The perturbations δij and δjk that apply to incoming and outgoing
arcs to and from node j are positively correlated, with a correlation coefficient that is
drawn uniformly in [0, 1].
As in the previous experiments, we solve each instance for various values
ν ∈ {0.15, 0.25, 0.35, 0.55, 0.75, 0.85, 0.95}. We evaluate each solution on a dataset that
contains 10,000 data points and compute the relative difference of objective value ρ as in
the previous experiments. Following the results obtained in the two previous experiments,
we fix the value K/m to 3 since it seems to result in a satisfactory trade-off between
accuracy and computational performance.
On all instances considered in these experiments, both models based on USVC and UAC
return the same solution for a given input parameter ν. Figure 3.12 shows the percentage
of optimal solutions found across all the instances with USVC and UAC . In the case of
USVC -based robust counterparts, about 55% of the instances are solved optimally within
the time limit of 10 hours. The final gap achieved on the remaining instances varies
from 0.42% to 34.22%, with an average value of 9.96%. On the other hand, all instances
are solved within 2 hours when the robust counterpart is defined using UAC , with 50%
and 75% reaching the optimal solution within 2 and 7 minutes, respectively. It is worth
noting that over the subset of instances solved to optimality with both uncertainty sets,
the solutions found were always identical.

5 Conclusion
In this paper we presented a method to build approximations of the SVC-based un-
certainty sets introduced in Shang et al. (2017) that preserves their polyhedral structure
while reducing significantly the number of additional variables and constraints in the re-
sulting robust counterpart. Our numerical experiments attest that the obtained solutions
with our approximate sets are comparable but require much less computational time (from
less than 1% to less than 50% in a few isolated instances) than the ones derived from the
original SVC-based uncertainty sets. A direct implication is that our procedure extends
the range of possible applications of the latter to problems of greater sizes. In addition,
any optimization procedure that relies on the polyhedral structure of the uncertainty set
may benefit from this approximation. From the practitioner standpoint, the parameter K

62
 Figure 3.12 – Proportion of TSP instances solved for USVC and UAC K/m = 3

used in the first phase of our method defines a trade-off between the quality of the approx-
imation and the reduction of the computational time. We test the influence of K on these
two aspects and provide some insights on appropriate values for this parameter. While
tuning its value may improve the performance in some specific applications, we observe
that K ≈ 3m, where m is the number of uncertain parameters, result in a good trade-off
in general. Finally, we notice that our approximation method is almost insensitive to the
type of uncertainty distribution in our experiments and to the value of the parameter ν
of the SVC algorithm. While our approximation does not provide the in-sample coverage
guarantee of the SVC-based uncertainty set (i.e. USVC contains at least a portion (1 − ν)
of D), our results show that for a sufficiently large value of K, the out-of-sample coverage
of both uncertainty sets are similar and the approximation does not significantly impact
the level of conservatism of the robust solution.
Despite these promising numerical results, we note that the use of a K-medoid pro-
cedure in its first phase makes our approach heuristic. Alternatively, one could think of
jointly optimizing the position and weight coefficient of the representatives of the support
vectors in a single phase. This results into large Mixed-Integer Non-Linear Programs that
are essentially more challenging and require new ideas to compute solutions efficiently.
Finally, notice that we rely on data points to approximate the original polyhedron ob-

63
tained by SVC. It will be interesting to study whether similar ideas are applicable to other
data-driven uncertainty sets such as the ones developed in Bertsimas et al. (2017) or Ning
and You (2018b).

64
 Chapter 4

A DATA -D RIVEN R OBUST O PTIMIZATION

APPROACH TO SUPPLY PLANNING

1 Introduction
In this chapter, we study a problem derived from a practical case encountered in the
aircraft industry. We consider the case of a single assembly line which manufactures
several products in order to satisfy the needs that are planned over a discrete and finite
planning horizon. Each final product has its own Bill Of Material (BOM) that involves
different components. The components are stored in a remote warehouse and they are
picked and delivered to the production line by a Third-party Logistic Provider (TPL)
following the manufacturer’s orders. In a sense, this problem is closely related to the
class of assemble-to-order systems, which are often encountered in the inventory control
literature. The need to define both component replenishment and production policies in
the presence of uncertain factors makes assemble-to-order systems particularly challeng-
ing. Large surveys on this topic have been presented in Song and Zipkin (2003) and more
recently in Atan et al. (2017).
In this chapter, we investigate the planning of picking orders when the available in-
formation related to picking times at the TPL is uncertain in a context where the time
available for these picking operations in a given period is bounded by a capacity. We model
this problem as an extension of the Multi-Item Capacitated Lot Sizing Problem (CLSP)
and propose a robust optimization approach to integrate the picking time uncertainties.
These picking time uncertainties are modeled as setup time uncertainties.
The capacitated lot sizing model is well-suited for production planning problems in
which (at least) one of the resources used in the production process is available in limited
quantity, imposing a capacity bound on the volume of goods produced during a period. In
practice, this can correspond to a maximum quantity of raw material, a maximum quan-

65
tity of produced items or a time capacity as in our case. The CLSP has been extensively
studied in the literature because of its relationship to a wide variety of practical prob-
lems arising in the industry. Several papers have reviewed the existing works on CLSP
an the different extensions of these models (Karimi and Fatemi (2003); Gicquel et al.
(2008)). The first application of such models in an uncertain context appears in Brandi-
marte (2006) who considers uncertain demand and proposes a multi-stage mixed-integer
stochastic programming model. Since then, most of the existing literature on CLSP under
uncertainty have focused on uncertain demand, including several RO approaches. Some
extensions of this basic setting to more complex versions of CLSP exist, see e.g. Coniglio
et al. (2018) who introduces two robust models for a CLSP with storage deterioration
or Abdel-Aal (2019) for a problem with setup times that allows overtime decisions. The
numerical experiments conducted in these studies show that the computation time needed
to solve robust formulations of a CLSP increases significantly with the number of peri-
ods in the planning horizon. To the best of our knowledge however, all existing works
on robust CLSP consider a budget-based uncertainty set and no study has focused on
uncertain setup times yet.
The contributions of this chapter are the following:

1. We provide a MIP formulation of a practical variation of a CLSP arising in the

aircraft industry.

2. We derive two robust formulations of the problem based on polyhedral uncertainty

sets. The first one uses the classical budget-based approach of Bertsimas and Sim
(2003, 2004) and the second one relies on the SVC-based uncertainty set of Shang
et al. (2017), or the approximate version introduced in the previous chapter.

3. We evaluate and compare the performances of the three models on simulated data.

The remainder of this chapter is organized as follows: Section 2 presents the supply
planning problem that motivates this study, along with a deterministic mathematical
formulation. In Section 3 we derive two tractable robust formulations of the problem
based on the budget-based and the SVC-based uncertainty sets, respectively. Finally,
Section 4 presents the experimental protocol and the numerical results obtained with
both robust formulations, while Section 5 contains our concluding remarks and some
perspective for future research directions.

66
 Warehouse Assembly line
Final product I

Components J

Figure 4.1 – Diagram of the problem

2 Problem statement and deterministic formulation

We study the case of a single assembly line that combines different components into
several final products over a finite and discrete planning horizon of T periods. Let [I]
and [J] the set of products and components, respectively. In each period t = 1, . . . , T ,
the system faces a demand dit for product i ∈ [I] that is assembled on demand (i.e. it is
impossible to manufacture a product in advance, then store it to satisfy a future demand).
For all product i ∈ [I] and component j ∈ [J], let rij be the number of components of
type j required to produce one unit of product i.
All components are assumed to be available at all time in a remote warehouse managed
by a TPL provider, who is in charge of delivering the assembly line based on its orders
for components, as represented on Figure 4.1. For each picking operation for component
j ∈ [J], the quantity collected by an operator is limited to a maximum batch size of mj
units. Note that several batches of the same components can be scheduled in the same
period. The picking time for a particular batch of component j of size qj can be divided
into two main parts. First, we consider a fixed time pj that corresponds to the travel time
between the shipping point and the zone where components of type j are stored, which is
spent regardless of the quantity that is transported. In addition, we consider a per-unit
picking time τj . Finally for all period t ∈ [T ], the total picking time spent by the TPL
provider to collect all the components delivered to the assembly line in period t cannot
exceed a maximum work capacity Ct .
Any demand for product i that is not satisfied immediately is backlogged until the

67
corresponding product is assembled in a subsequent period. Each such unit incurs a
per-unit backlogging penalty cost bi for each period of delay. Whenever a component j is
available on the assembly line but is not immediately used to manufacture an end product,
it disturbs the production process by interfering with people and other goods moving
nearby. We model this situation with a per-unit, per-period obstruction cost oj . The
problem consists in planning the quantity of each component delivered to the assembly
line by the TPL in each period such that the sum of the obstruction and backlogging
costs is minimized. We let qjt ∈ R+ be the quantity of components j ∈ [J] that is ordered
by the assembly line for a delivery in period t. These quantities induce the number of
distinct picking operations of components j ∈ [J] performed during period t represented
as the integer variable xjt ∈ N. For any t ∈ [T ], the variable uit denotes the number of
final product i ∈ [I] produced during period t and Ot corresponds to the total obstruction
cost for period t. At the end of each period t ∈ [T ] the stock level sjt of each component
j ∈ [J] is given by:
t
X X
sjt = sj0 + (qjt − rij uik ) (4.1)
k=1 i∈[I]

We formulate the problem with the following MIP:

 
T
X X t
X
min Ot + bi (dik − uik ) (4.2)
t=1 i∈[I] k=1
 
X t
X X
s.t. Ot ≥ oj sj0 + (qjt − rij uik ) ∀t ∈ [T ] (4.3)
j∈[J] k=1 i∈[I]
t
X t
X
uik ≤ dik ∀i ∈ [I] , ∀t ∈ [T ] (4.4)
k=1 k=1
X
(pj xjt + τj qjt ) ≤ Ct ∀t ∈ [T ] (4.5)
j∈[J]

xjt ≥ ⌈qjt /mj ⌉ ∀j ∈ [J] , ∀t ∈ [T ] (4.6)

qjt , Ot ∈ R+ ∀j, ∈ [J] , ∀t ∈ [T ] (4.7)
uit , xjt ∈ N ∀i ∈ [I] , ∀j ∈ [J] , ∀t ∈ [T ] (4.8)

The objective (4.2) aims at minimizing the total cost incurred over the planning hori-
zon. Note that the backlogging term is always positive due to constraints (4.4), which
ensure that the production plan up to a given period cannot anticipate demands in future
periods. Constraints (4.3) compute the total obstruction cost for period t ∈ [T ]. Con-

68
straints (4.5) ensure that the capacity of the logistic provider is respected by the picking
plan during a given period t ∈ [T ]. Finally, constraints (4.6) link the quantity of compo-
nent ordered qjt and the number of batches xjt , while constraints (4.7) and (4.8) define
the domain of the decision variables.

3 Robust models for uncertain setup times

In practice, the manufacturer often has either incomplete or imprecise knowledge on
the picking times. As a consequence, some combinations of his orders may exceed the
picking capacity of the TPL provider, forcing the latter to postpone some operations
to subsequent periods. This delay in the delivery of some components is susceptible to
induce two types of inefficiencies on the assembly line as missing components (i) prevent
the manufacturer to assemble some of the products, leading to backlogging penalty costs
and (ii) leave the other components in the BOM on the border of the line, disturbing
other production operations. This double negative effect strongly incentivizes the planner
to protect his decisions against uncertain setup times. In what follows, we assume that
we have access to a set of historical setup times P = {p(1) , . . . , p(N ) } in order to tune the
different parameters of the uncertainty sets.
As the uncertain parameters of this problem are the setup picking time p, its robust
formulation only affects the capacity constraints (4.5). Given uncertainty set Ut for a
period t ∈ [T ], the robust counterpart of this constraint is expressed as follows:
X
pj xjt + τj qjt ≤ Ct ∀p ∈ Ut , ∀t ∈ [T ] (4.9)
j∈[J]

which is equivalent to the following constraint (see Chapter 1):

 
X  X
max pj xjt + τj qjt ≤ Ct , ∀t ∈ [T ] (4.10)
p∈Ut  
j∈[J] j∈[J]

3.1 Budget-based robust formulation

As a first robust approach, we consider the classical uncertainty set of Bertsimas and
Sim (2003, 2004) already defined in previous chapter and defined by :

UtΓ = {p |pjt = p̄j + p̂j zjt | ∥z∥1 ≤ Γt , ∀j ∈ [J] } (4.11)

69
, where p̄j is the average value of pj evaluated on P and p̂j is a maximal deviation. The
scaled deviation coefficient zjt = (pj − p̄j )/p̂j ∈ [−1, 1] ensures that pjt ∈ [p̄j − p̂j , p̄j + p̂j ]
for all components j ∈ [J] and period t ∈ [T ] and that the total scaled deviation of p
never exceeds the budgets of uncertainty Γ = (Γ1 , . . . , ΓT ).

Proposition 1. The robust counterpart of problem (4.2)-(4.8) based on the uncertainty

sets UtΓ , t ∈ [T ] is equivalent to the following MIP, after introducing the additional
J×T
variables v ∈ RT+ and w ∈ R+ :

 
T
X X t
X
min Ot + bi dik − uik  (4.12)
t=1 i∈[I] k=1

s.t. (4.3), (4.4), (4.6) − (4.8) (4.13)

X X
p̄j xjt + τj qjt + vt Γt + wjt ≤ Ct ∀t ∈ [T ] (4.14)
j∈[J] j∈[J]

p̂j ≤ vt + wjt ∀j ∈ [J] , ∀t ∈ [T ] (4.15)

vt , wjt ∈ R+ ∀j ∈ [J] , ∀t ∈ [T ] (4.16)

Proof. We formulate the robust capacity constraint (4.10) as:

 
X X  X
p̄j xjt + max p̂j xjt zjt + τj qjt ≤ Ct , ∀t ∈ [T ] (4.17)
j∈[J] zt ∈[0,1]|J| 
j∈[J]

j∈[J]
∥zt ∥1 ≤Γt

Considering an optimal solution x∗ to the problem, for each t ∈ [T ], the inner maximiza-
tion problem is equivalent to a linear optimization problem. As values in x are all positive
or null, we can restrict zt ∈ [0, 1]|J| and thus replace the l1 -norm ∥zt ∥1 by a simple sum
P
j∈[J] zjt . This gives us the following equivalent LP:

p̂j x∗jt zjt

X
max
j∈[J]
X
s.t. zjt ≤ Γt (4.18)
j∈[J]

0 ≤ zjt ≤ 1, ∀j ∈ [J] .

After introducing dual variables vt and wjt for j ∈ [J], we obtain the dual formulation of

70
the problem:
n
X
min vt Γt + wjt
j∈[J]

s.t. vt + wjt ≥ p̂j x∗jt ∀j ∈ [J] (4.19)

vt ≥ 0
wij ≥ 0 ∀j ∈ [J]

Since the primal problem (4.18) is feasible and bounded, its dual problem (4.19) is also
feasible and bounded and their optimal values coincide. Therefore, we can replace the
inner maximisation problem of period t ∈ [T ] by the objective function of (4.19), with the
associated constraints. This gives us the (linear) robust counterpart (4.12)-(4.16).

3.2 SVC-based robust formulation

As a second robust formulation, we propose to use the SVC-based uncertainty set of

Shang et al. (2017) presented in the previous chapter. Let us recall that the SVC approach
applied on the set of historical picking times P define a subset S of data points, called
support vectors (SV), that lie on or outside of the boundary of the uncertainty set. A
subset B ⊆ S of these points, referred to as boundary support vectors (BSV), are located
exactly on the boundary of the uncertainty set. The definition of this uncertainty set
ensures that B ⊆ S ⊆ D and |S \ B| ≤ N ν ≤ |S|. In addition, the SVC algorithm
also computes for each historical data point l ∈ [|D|] a value αl ∈ [0, 1/νN ] such that
αl = 0 iff l ∈ D \ S, αl = 1/νN for all l ∈ S \ B and 0 < αl < 1/νN when l ∈ B. Let
α = [α1 , . . . , αN ] be the vector of obtained coefficients, Shang et al. (2017) use α to define
the following data-driven uncertainty set:
 
 
UtSV C = p αl ∥Q(p − p(l) )∥1 ≤ θ 
X
(4.20)
l∈S

1
P ′

where Q = Σ− 2 whith Σ the covariance matrix of p and θ = min
′ l∈S αl ∥Q(u(l ) − u(l) )∥1 .
l ∈B
|J|×|S|
Proposition 2. Using for all ∀t ∈ [T ] the additional variables λt , µt ∈ R+ and
ηt ∈ R+ leads to the following MIP formulation for the robust counterpart of problem
(4.2)-(4.8) based on the uncertainty sets UtSV C :

71
  
X X t
X
min Ot + bi (dik − uik ) (4.21)
t∈[T ] i∈[I] k=1

s.t. (4.3) − (4.8) (4.22)

(µtl − λtl )T Qu(l) + ηt θ +
X X
τj qjt ≤ Ct ∀t ∈ [T ] (4.23)
l∈S j∈[J]
X
Q(λtl − µtl ) + xt = 0 ∀t ∈ [T ] (4.24)
l∈S

λtl + µtl = ηt αl 1 ∀l ∈ S, ∀t ∈ [T ] (4.25)

λtlj , µtlj , ηt ∈ R+ ∀j ∈ [J] , ∀l ∈ S, ∀t ∈ [T ] (4.26)

Proof. Considering an optimal solution x∗ to the problem, the inner maximization prob-
lem of constraint (4.10) for a given t ∈ [T ] is equivalent to the following linear optimization
problem:

max pT xt
αl 1T vl ≤ θ
X
s.t. (4.27)
l∈S

− vl ≤ Q(p − p(l) )≤ vl ∀l ∈ S

|J|×|S|
We introduce the dual variable λt , µt ∈ R+ and ηt ∈ R+ , ∀t ∈ [T ] to define the dual
of the problem:

(µtl − λtl )T Qu(l) + ηt θ

X
min
l∈S
X
s.t. Q(λtl − µtl ) + xt = 0 ∀t ∈ [T ] (4.28)
l∈S

λtlj , µtlj , ηt ∈ R+ ∀j, ∈ [J] , ∀l ∈ S, ∀t ∈ [T ]

Since the primal problem (4.27) is feasible and bounded, its dual problem (4.28) is also
feasible and bounded and their optimal values coincide. Therefore, we can replace the
inner maximisation problem of period t ∈ [T ] by the objective function of (4.28), with
the associated constraints. The (linear) robust counterpart (4.21)-(4.25) follows.

72
4 Experimental results
The aim of this section is to evaluate the performances of the proposed optimization
models, both in terms of their practical computational tractability and the quality of the
solutions obtained. To do this, we generated a set of instances with various configuration
(i.e. number of products, number of components, number of periods). In what follows,
we will use the notations Det to denote the deterministic models or their solutions, U Γ to
denote the budget based robust models or their solutions and U SV C and U AC to denote the
SVC-based models and its approximation introduced in the previous chapter, respectively.
All models are built from historical data sets P of 1000 data points.

Instance generation We first describe several rules that we apply to generate the
instances used to evaluate the performances of the different models presented above. The
specific assumptions we consider are listed below, with the objective to reflect the original
use case from the aircraft industry that motivates this work: rij is chosen randomly from
an interval [rj , rj ], the smallest interval is [2,4] and the widest is [10,25]. Each product i
requires a subset Ji ⊆ [J] of components. The obstruction cost oj of component j is chosen
P
randomly in [0, 1], while bi = ρ j∈[J] rij oj is a multiple of the sum of the obstruction cost
of its components. Here ρ represents a ratio between the penalty incurred for backlogging
a demand for product i versus the maximum obstruction cost of its components. In our
numerical experiments, we use ρ = 5, but this value can be adjusted depending on the
context. In order to primarily focus on the impact of the setup times uncertainty, we
disregard the linear part of picking times and set τj = 0 and mj = 3 maxi∈[I] rij for all
j ∈ [J]. Demands are deterministic with values drawn from a uniform distribution whose
range vary among products. We compute the capacity of the TPL Ct from the average
time needed to pick the required components in each period, that we increase by a given
percentage to ensure that it is not always saturated.

Setup time distribution We assume that the TPL provider organizes the storage of
the components in order to optimize the picking operations. Specifically, components are
stored in such a way that their accessibility improves with their order frequency. We
consider three types of components and separate them based on their storage area. We
assume that both the mean and the variability of the setup times decrease with component
accessibility. In our instances, we thus generate setup times using Gamma distributions
with three different shape parameters, one for each type of component. We also assume

73
that picking times are influenced by the behaviour of the operators and can be correlated if
operators groups orders of components or use particular picking routes. Thus, we generate
random correlations among picking time of components in the same storage area.
In practice, setting the budget parameter Γt is hard and the performance obtained for
a given parameter value heavily depends on the characteristics of the instance. In order to
compare both models, we set Γt to cover a certain portion of P, by iteratively increasing
its value until we reach the desired coverage.

Evaluating solutions As we are considering uncertain parameters, we can not com-

pare solutions only based on the objective value of the optimization models directly. In
addition, we cannot rely on a simple evaluation of the constraint satisfaction as the main
objective of this problem is to increase the satisfaction of the final product demand with
a reasonable increase of the amount of components held on the side of the assembly line.
Thus in order to evaluate the quality of solutions obtained, we rely on a simulation process
following four steps:
1. Generate a set of 104 setup time scenarios S following the same distribution than
P. Each scenario is composed of one setup time for each component j ∈ [J] and
each period t = 1, . . . , T . This unique set will be used to evaluate all the solutions
for a given instance.
2. For a given solution x∗ , q ∗ simulate the picking operations at each period with each
scenario in S to obtain the real quantities of components delivered on the assembly
line. This simulation take into account the reported quantity if the capacity of a
given period is exceeded.
3. Compute the optimal production plan for the quantity delivered by the simulation.
4. Evaluate the resulting production plan, the costs, the satisfaction of the demand
and the amount of components held next to the assembly line.
As in Chapter 3, all models are implemented in Julia 1.9 with JuMP and solved with
CPLEX 20.1 using a single thread on a Linux, Ubuntu 20.04.2 LTS, equipped with an
Intel Xeon Gold 6230 clocked at 2.10GHz.

4.1 Scalability of the SVC-based model

Recall from the previous chapter that the number of variables and constraints of
the SVC-based uncertainty set is proportional to the number of uncertain parameters of

74
 the problem and the number of data point in the training sample. This is particularly
problematic in the problem considered in this chapter, which includes multiple constraints
involving uncertain parameters. In this first set of experiments we consider instances with
2 final products, 5 components and a planning horizon of T = 7 periods. Table 4.1 shows
the average and 0.95 quantile of solution costs, as well as the average computation time
of both models when the uncertainty sets are tuned to cover different portion of P. As we
consider a set of small instances, we restricted the computation time to 1 hour for both
models and report the remaining CPLEX gap (in percent) when no optimal solution is
found within the time limit. In this table, we partition the results into three cases based
on the desired coverage of P:

Case 1 When the coverage is 0.15 or 0.35 both models are solved to optimality and pro-
vide comparable performances in term of costs. However, using the approximate
uncertainty set U AC reduces the computation time by more 99%.

Case 2 When the coverage is 0.55, the original SVC-based formulation presents an average
CPLEX optimality gap of 3.54% after 1 hour while our approximation converges to
an optimal solution within 10.47 seconds on average. The impact of this gap on the
solutions seems significant as the average cost is increased by more than 6% and the
0.95 quantile by more than 13% when using U SV C .

Case 3 When the coverage is high (0.75 or 0.95) both models are extremely conservative,
leading to solutions with high average costs that are nearly unaffected by the un-
certainty as the difference between the average value and the 0.95 quantile suggests.
Moreover, this comes with a significant increase of the CPLEX CPU time that is
not influenced as much by the size of the models and affects U AC in a similar fashion
as U SV C .

We restricted the CPLEX CPU-time to 1 hour as our objective is to be able to obtain

good solutions in a reasonable amount of time. The results show that both models lead
to comparable solutions in term of costs but with a significant reduction of the CPLEX
CPU-time, in particular on small instances for which both U SV C and U AC find an optimal
solution on all instances. This suggests that U AC may be applied to larger instances to
obtain solutions in a reasonable amount of time. They illustrate that the approximate
uncertainty set designed in Chapter 3 gives an opportunity to the practitioners to use the
information they can extract from the data available through an unsupervised learning
procedure without the computational burden that usually comes with it.

75
Table 4.1 – Comparison of U SV C and the approximation U AC on small instances with 2
final products, 5 components and time horizon T = 7 periods.

U SV C U AC
Coverage
Average Q0.95 time (s) gap Average Q0.95 time (s) gap
0.15 1487 4050 543.93 0 1485 4015 1.65 0.0
0.35 1324 3309 1112.9 0 1319 3275 6.41 0.0
0.55 1329 2727 3600.88 3.54 1236 2362 10.47 0.0
0.75 3508 3842 3600.88 42.5 3164 3612 952.41 0.0
0.95 9451 9446 3291.47 0 9451 9446 2199.95 0.0

4.2 Comparing U Γ and U AC

In all the robust models we present, we use the portion of data coverage as a parameter
to define the degree of robustness of the solution. The impact of this parameter on the
cost of solutions is represented in Figure 4.2,which compares the evolution of the average
(4.2a) and 0.95-quantile (4.2b) costs obtained by both U Γ and U AC when the portion of
P covered increases. The instance includes 5 final products, 10 components and T = 14
periods. For both models we observe that increasing the portion of historical data covered
significantly reduces the total cost in extreme cases. Figure 4.2b shows that when the
coverage increases, the 0.95-quantile of costs is reduced from more than 25000 to less
than 17500 for U Γ and to less than 14000 for U AC . On the other hand, this improvement
for extreme cases comes with a deterioration of the average solution cost, as shown in
Figure 4.2a. In the case of U Γ the average cost decreases from 10600 to 9900 as the
coverage grows up to 0.25, but increases with the coverage afterward. While we observe
the same behavior when using U AC instead, the decrease in average cost (from 9900 to
8000) continues until a coverage value of 0.35, but starts going up beyond that point. For
all coverage values tested, U AC clearly outperforms U Γ by providing lower average and
extreme costs. While the two 0.95-quantile average costs are comparable for coverage
values between 0.2 and 0.35, this comes at the expense of a steep raise in average cost
for U Γ . We observe a similar phenomenon affecting U AC , but starting at higher values of
coverage and with milder increasing slope of average costs.
Table 4.2 details the obstruction and backlogging costs obtained by both methods
on a set of 10 instances with 5 final products, 10 components and a planning horizon
of T = 7 periods when the coverage of P varies from 0.1 to 0.3. The first observation

76
 (a) Average cost (b) 0.95-quantile cost

Figure 4.2 – Evolution of the average cost and 0.95-quantile obtained with U Γ and U AC
on an instance with 5 components 10 products and planning horizon T = 14 periods.

is that U AC simultaneously provides lower obstruction and backlogging costs when the
coverage is 0.1 and 0.2, which explains the difference observed on Figure 4.2. When the
coverage is 0.3, U Γ achieves a greater reduction of the obstruction cost at the expense
of a significant increase in backlogging costs, which suggests that its representation of
possible parameter values includes unlikely scenarios which do not lead to anticipate
enough the need for some components to satisfy the demand on time. In comparison,
U AC provides a more significant reduction of backlogging costs with higher obstruction
costs and lead to lower costs overall. In the last line of 4.2, when the coverage is 0.5, both
approaches provides over-conservative solution increasing simultaneously the obstruction
and backlogging costs. This conservatism is emphasized by the static decision context
and may be reduced in a more dynamic one where the decisions can be adjusted to
the realization of uncertain parameters over time. Table 4.3 shows the evolution of the
average CPLEX CPU-time and the CPLEX optimality gap after 1 hour for the same set of
instances as the coverage of P increases. We observe that for both models, increasing the
coverage increases the CPLEX CPU-time, for example the solutions of U AC are obtained
in 2 seconds on average when the coverage is 0.1 and in 193 seconds when the coverage
is 0.3. In addition, when the coverage is high (≥ 0.5) the CPLEX CPU-time increase
significantly for both models (2200 seconds for U Γ ) and the optimal solution is never
obtained in 1 hour for U AC .
Table 4.4 compare the obstruction and backlogging costs obtained by both methods
on a set of 10 instances with 10 final products and 20 components for different coverage

77
Table 4.2 – Comparison of U Γ and the approximation U AC on a set of 10 instances with
5 final products, 10 components and time horizon T = 7 periods.

UΓ U AC
Coverage Obstruction Backlogging Obstruction Backlogging
Average Q0.95 Average Q0.95 Average Q0.95 Average Q0.95
0.1 7022 8137 4286 16852 −10.89% −13.25% −17, 49% −10.72%
0.2 8677 9616 3251 14600 −17.76% −17.73% −6.21% −9.76%
0.3 7985 8654 6064 12502 +34.41% +30.03 −71.82 −34.32
0.5 9125 10189 13116 13782 −10.17% −11.31 −24.41 −22.27

Table 4.3 – Average CPLEX CPU-time and optimality gap of U Γ and the approximation
U AC on a set of 10 instances with 5 final products, 10 components and time horizon T = 7
periods.

UΓ U AC
Coverage
time (s) gap (%) time (s) gap (%)
0.1 1 0 2 0
0.2 3 0 12 0
0.3 9 0 193 0
0.4 12 0 675 0
0.5 2200 2 3600 5
0.6 3600 7 3600 9

of P. In line with precedent results of Table 4.2, this table show that U AC outperform
U Γ by providing solutions with lower average and extreme costs. However, the relative
difference seems to be smaller for this set of larger instances.

5 Conclusion
In this chapter, we consider a SVC-based robust optimization approach to plan the
picking operations of an assembly line when the latter are subcontracted to an outside
service provider. Through a set of experiment we showed that using the classic SVC-
based approach of Shang et al. (2017) is not a viable option to solve this problem in
practice due to its long solving time. We propose to instead rely on the approximation
introduced in Chapter 3 to obtain solutions quickly. We compare these solutions with
a classical budget-based approach similar to the one introduced in Bertsimas and Sim
(2003, 2004). Experimental results show that the proposed method efficiently reduces the

78
Table 4.4 – Comparison of U Γ and the approximation U AC on a set of 10 instances with
10 final products, 20 components and time horizon T = 7 periods.

UΓ U AC
Coverage Obstruction Backlogging Obstruction Backlogging
Average Q0.95 Average Q0.95 Average Q0.95 Average Q0.95
0.1 14343 17284 9674 25354 −12.28% −14.17% −13, 47% −9.58%
0.2 17192 20017 7191 22405 −9.51% −12.18% −9.21% −10.40%
0.3 15030 18195 13011 20610 +17.42% +15.78% −52.15% −24.89%

impact of uncertainties on the performance of the assembly line. By controlling the level
of robustness of the model, good tradeoffs can be found between reducing the backlog
of assembly operations and minimizing the quantity of components stored on the border
of the assembly line. In comparison with a classical budget-based robust optimization
model, the data-driven model yields better performance on test instances with up to 10
different final products and 20 different components. The comparisons between solutions
provided by both methods shows that the data-driven approach can lead simultaneously
to a higher satisfaction of the demand (smaller backlogging cost) and a lower level of
obstruction of the assembly line.
Despite the computational improvement provided by the use of the approximation,
using an SVC-based approach lead to a significant increase of the CPLEX CPU-time
compared to the budget based method. In addition, the computation time needed to
solve both types of models increases significantly with the desired level of robustness
(i.e coverage of historical data). An interesting perspective to apply this approach on
larger industrial instances would be to combine the proposed approximate SVC-based set
with adapted heuristic approaches, e.g. the well-known Relax-and-Fix, Fix-and-Optimize
one, to provide good quality solutions in a reasonable amount of time. Notice that since
robust models are reformulated as a deterministic robust counterpart, existing heuristic
approaches from the literature on CLSP could be applied without further modifications.
An other relevant perspective for this specific application would be to consider a
dynamic decision context in which the quantity ordered to the supplier can be adjusted
dynamically to take the reported picking operations into account. This would lead to a
better estimation of the consequences of picking time uncertainties on the assembly line
and would reduce the cost of robust solution. However, the decentralized decision-making
framework we consider assumes that picking decisions are made by the logistic provider,

79
which makes the description of the problem dynamics particularly challenging. Whether
the TPL policy can be included in the uncertainty set or incentivized through the orders
placed by assembly line manager remain an open question that should be addressed if one
wishes to derive a formulation that faithfully relates to practical situations.

80
 Chapter 5

DATA -D RIVEN R OBUST INVENTORY

MANAGEMENT UNDER UNCERTAIN
DEMAND AND LEAD TIME

1 Introduction
In the two previous Chapters, we studied the efficiency of a non-parametric DDRO
approach (SVC-based) that can be used as a black-box on any type of uncertainty without
further tuning effort when data are available. In this Chapter, we propose to formulate our
uncertainty sets from historical data but we assume that random parameters are driven
by a given stochastic process (ARIMA) or have specific relationships with each others
(Linear correlations). Moreover, while the SVC-based uncertainty set was used to model
a single type of uncertainty in the previous models (item weights, demands, travel times,
picking times) in this chapter, we study a more complex problem in which two types of
uncertain parameters (demand and lead time) are modeled simultaneously.
Uncertain demand and lead time are two of the main factors driving up inventory man-
agement costs. Poor estimation of demand can lead to overstocking, generating additional
costs, or to stock-outs that can have negative effects on future sales and profitability. In
the same way, extended lead times or delivery delays increase the risk of stock-outs if
the inventory level is insufficient. One may choose to adopt a conservative approach in
this case and overestimate lead times but this can significantly increase storage costs. It
has been pointed out in the inventory management literature that the majority of the
existing works focus on demand uncertainty while uncertainties on delivery or produc-
tion lead times are seldom considered Moreover, most contributions focus on one type of
uncertainty even if they are likely to occur simultaneously in practical problems (Ben Am-
mar et al., 2013). The majority of research papers assume that the true distribution of

81
uncertain parameters is either known or well approximated by fitting the first empirical
moments of the sample available. They then compute optimal reorder points and order
quantities based on this information. However, this situation is often unrealistic as the
true distribution is generally unknown and hard to characterize in practical cases. In
addition, errors in the evaluation of this assume distribution of lead time can lead to
sub-optimal ordering policies (see e.g. Eppen and Martin (1988) who show that such
approaches are sub-optimal in specific cases).
To overcome this drawback, Scarf’s method paved the way for distributionally robust
optimization with application to inventory management (Scarf, 1958; Gallego, 1992; Moon
and Gallego, 1994). More broadly, RO has emerged in recent years as an appropriate
approach to address various inventory management problems, especially in a distribution-
free setting (see Chapter 1 for a short review of existing references on the topic).
It is a common practice for inventory or production managers to adjust their replen-
ishment or planning decisions based on demand forecasts, often obtained through time
series models that are tuned on historical data to faithfully capture the dynamics of a
demand process. The classical methodology applied both by researchers and practitioners
aims at designing a model that produces the most accurate predictions before using them
as input parameters for their decision support tools. Among the options available, Auto-
Regressive Moving Average (ARMA), or its integrated version Auto-Regressive Integrated
Moving Average (ARIMA), process is a popular one to produce high-quality forecasts
(Udom and Phumchusri (2014); Syntetos et al. (2008)). The natural relationship that
link better predictions accuracy and improved decisions is not always true, as pointed
out by several studies. Most of the works that apply robust optimization to inventory
management and lot-sizing problems make the standard (but restrictive) assumption of in-
dependence of uncertain parameters. However, multiple papers have shown that demand
processes may display autocorrelation (Erkip et al., 1990; Disney et al., 2006), in some
cases with significant impact on inventory policies (Luong, 2007; Charnes et al., 1995).
In addition, existing contributions on supply chain management have stated that lead
times may be correlated with the workload of suppliers (Pahl et al., 2005, 2007). Indeed,
a succession of orders may saturate a supplier’s production capacity and thus increase
the delivery lead time and variability. This correlation could also be strengthened in a
competitive context if several retailers are facing high demands in the same period. On
the contrary, a sequence of periods with low demand may incentivize suppliers to group
production or transportation of items, inducing an increase in lead time. One way to

82
reduce the negative impact of load dependent lead times on the inventory management
problem is to consider a joint production and inventory problem (Boute et al., 2014) or to
use the workload of the supplier as an input to estimate the delivery lead time. However,
in practice, this information is not always available for retailers and the ordering policy
is often based on the observed delivery lead times. In this Chapter, we consider a robust
multi-period, single-item inventory management problem with uncertain demand and lead
time where all ordering decisions should be made before the beginning of the time hori-
zon (Static Uncertainty model as described by Bookbinder and Tan (1998)). To address
this problem, we provide an easy to implement, distribution-free, methodology that relies
on existing tools and reduces the conservatism of classic robust optimization approaches
when the uncertain parameters are not identically and independently distributed (i.i.d).
In addition, we extend (Thorsen and Yao, 2017) in different directions: 1. We assume that
the demand is correlated over time and propose a Data Driven Uncertainty Set based on
autoregressive processes to capture this correlation. 2. We assume that the demand is
correlated with the delivery lead time and propose a joint uncertainty set for demand and
lead time to handle both types of uncertainty. 3. We consider fixed ordering costs and a
discrete lead time uncertainty set. We conduct a set of experiments to demonstrate that
our data-driven uncertainty sets outperforms classic uncertainty sets by integrating the
information retrieved from the available historical data.
The rest of the Chapter is divided as follows: Section 2 presents the mathematical
model and its robust formulation. Section 3 presents the proposed Data-Driven uncer-
tainty sets. Section 4 presents the computational experiments settings, the computation
results and discussions.

2 Problem and robust solution approach

We study an inventory Lot-Sizing problem over a discrete and finite horizon of T
time periods. The goal is to minimize the costs incurred to manage the inventory of a
single item in order to satisfy an uncertain sequence of demands (d1 , . . . , dT ) that follows
a specific (but unknown) stochastic process. At the beginning of each period t ∈ [T ],
the manager may place an order for a quantity of qt ≥ 0 units, incurring an ordering
cost of s1{qt ≥0} + cqt , where s is a fixed setup cost and c is a per-unit ordering cost. We
consider that an order does not replenish the inventory immediately and instead assume
that the units ordered are received after a discrete lead time of lt ∈ [1, L] periods. When

83
the available inventory in period t is insufficient to cover demand dt , the unmet demand
is backlogged to be satisfied in a subsequent period. At the end of each time period, the
system also incurs a holding cost h for every unit physically held in the inventory and a
backlogging penalty cost b for every unit of postponed demand.

2.1 Deterministic model

In a deterministic model, both lead time and demand vectors are known. For all k ≤ t,
we model lead times using binary parameters

1

if the quantity qk ordered in k is delivered in or before period t
δkt =  (5.1)
 0 otherwise

The inventory level at the end of period t is given by the following equation:

t
X
It = I0 + (δkt qk − dk ) (5.2)
k=1

where I0 is the initial inventory level. The corresponding inventory cost incurred in period
t is thus yt = max{hIt , −bIt }.
We use the above observation to formulate the problem as the Mixed Integer Program
(MIP) below :

T
X
min sxt + cqt + yt (5.3)
x∈{0,1}T t=1
q,y∈RT
+
t
!
X
s.t. yt ≥ h I0 + (δit qi − di ) ∀t ∈ [T ] (5.4)
i=1
t
!
X
yt ≥ −b I0 + (δit qi − di ) ∀t ∈ [T ] (5.5)
i=1

qt ≤ M x t , ∀t ∈ [T ] (5.6)

where x ∈ {0, 1}T are binary variable such that xt = 1 if an order is placed in period t
and 0 otherwise and M is a sufficiently large value defined such that constraint (5.6) is
P
always true if xt = 1, e.g. M = t∈[T ] dt .

84
2.2 Robust optimization model

Following the classical paradigm of robust optimization, we consider that both the
demand vector d and delivery matrix δ belong to uncertainty sets Ud and Uδ , respectively,
similar to the ones introduced by Thorsen and Yao (2017) for an inventory problem with
no fixed costs. The goal is then to optimize against the worst parameters values in these
sets by solving the following MIP:

T
X
min T max minT sxt + cqt + yt (5.7)
x∈{0,1} d∈Ud y∈R
δ∈Uδ
+ t=1
q∈RT+
t
!
X
s.t. yt ≥ h I0 + (δit qi − di ) ∀t ∈ [T ] (5.8)
i=1
t
!
X
yt ≥ −b I0 + (δit qi − di ) ∀t ∈ [T ] (5.9)
i=1

qt ≤ M x t ∀t ∈ [T ] (5.10)

The objective function (5.7) minimizes the maximum value of the sum of fixed/unit or-
dering costs and holding/backlogging costs over the planning horizon. Constraints (5.8)
and (5.9) define the holding/backlogging cost yt due to storage or stockouts at the end
of period t ∈ [T ]. Finally, constraints (5.10) link the binary vector x with the vector of
ordered quantity q.
As mentioned in Chapter 1, a popular choice in robust inventory management is the
budget-based uncertainty set introduced by Bertsimas and Sim (2003, 2004). It has been
used in Bertsimas and Thiele (2006) to define admissible values for a demand vector from
so-called budgets of uncertainty Γ:

t
( )
UdΓ = d = d¯ + dˆ⊺ z m
X
zi ≤ Γt ∀t ∈ [T ] , z ∈ [−1, 1] (5.11)
i=1

where d¯ is the vector of nominal demand and dˆ = [dˆ1 , . . . , dˆT ] is the vector of maximum
deviation from the nominal demand. The tuning parameters Γt control the degree of
robustness of the model by bounding the cumulative deviation of the demands from their
nominal values. The value of d¯ and dˆ define how UdΓ covers the possible values of d, which
may have a significant influence on the solutions obtained. Hence they must be selected
carefully to include a predefined fraction of the demands, e.g. 95% (see Bertsimas and

85
Thiele (2006)). More recently, Thorsen and Yao (2017) and Thevenin et al. (2022) pro-
posed a robust approach to inventory management problems under lead time uncertainty
by introducing an additional uncertainty set on δ:
 


 δtk ≤ δtk+1 , ∀t = 1, . . . , T − 1, t ≤ k ≤ T 



 



 T X
X T 


(1 − δtk ) ≤ Γ

 

 
UδΓ = δ t=1 k=t (5.12)

 

δtt+L = 1, ∀t = 1, . . . , T

 


 


 

 
δtk ∈ [0, 1], ∀t = 1, . . . , T, t ≤ k ≤ T

 


where L is an upper bound on the lead time value and the budget of uncertainty Γ bounds
the cumulative lead time over the planning horizon. Clearly δ must be a matrix of binary
values to ensure that lt = Tk=t (1 − δtk ) is an integer value for all t = 1, . . . T . Note
P
2
that Uδ relaxes this constraint by considering instead that δ ∈ [0, 1]T to obtain a convex
uncertainty set. This models a continuous lead time, effectively allowing an order to be
delivered in more than one single period. For instance, a vector δ1 = [0.5, 0.5, 0.5, 1, 1]
corresponds to a scenario in which half of the order q1 is received in period 1 and the
second half is received in period 4.

3 Data-driven uncertainty set

In this section, we apply time series analysis on available data to take advantage of
the structure of classical stochastic processes used to model demands. In particular, we
exploit properties of the type of processes that is often used as a demand forecasting tool
to incorporate correlation information between parameters and propose an alternative
uncertainty set that is less conservative than UdΓ . This approach differs from previous
data-driven uncertainty sets such as the SVC based uncertainty set studied in Chapter 3
and 4 as we assume that the demand follows a specific stochastic process. After presenting
the theoretical foundations to build this uncertainty set, we provide insights on fitting the
different parameters from historical data and propose an extension that incorporates the
influence of past demands on order lead times (production and/or transportation).

86
3.1 Data driven uncertainty set based on ARMA(p,q) process

Many existing works on inventory management rely on time series analysis to produce
forecasting models for demand process (see e.g. Duc et al. (2008) and included references).
Among the most popular ones, Autoregressive Moving Average (ARM A) processes often
offer enough flexibility to derive good quality predictions, with the objective to make
better informed production and/or ordering decisions. Besides its linear nature that
facilitates integration into MIP formulations, another key advantage is that one can easily
fit demand data to such processes using one of the many existing off-the-shelf libraries
dedicated to time series analysis.

3.1.1 General ARM A(p, q)-based uncertainty set

The demand (Dt )t∈N is an ARM A(p, q) process of mean 0 if it satisfies the following
relationship:

Dt = 0 ∀t ≤ 0
Dt + β1 Dt−1 + · · · + βp Dt−p = ϵt + θ1 ϵt−1 + · · · + θq ϵt−q otherwise

where β ∈ Rp and θ ∈ Rq are the vectors of autoregressive and moving-average coefficients

and where ϵ is a vector of i.i.d. noise of mean 0 and variance σϵ2 . If the process has mean
µ ̸= 0, the process defined for all t as Dt − µ satisfies the above equations.
Given a vector of N + 1 consecutive (historical) demand values [d−N , . . . , d0 ] that are
assumed to be generated by an ARM A(p, q) process, one can find the Maximum Like-
lihood Estimators (MLE) of parameters µ, σ, β1 , . . . , βp , θ1 , . . . , θq assuming a normally
distributed random noise ϵt .One can then use classical results on ARM A processes to ob-
tain the prediction interval for a given confidence level 1 − α. We define the ARM A(p, q)
based uncertainty set as:
 p
X q
X

dt = β0 + βk dt−k + θk ϵt−k + ϵt ∀t ∈ [T ] 

 

 

 
k=1 k=1

 


 

t √

 

 X 
d | ϵj | ≤ Γϵ σϵ t ∀t = 1, . . . , T
 
ARM A(p,q)
Ud = (5.13)

 j=1 


 

dt ≤ dt ≤ dt ∀t ∈ [T ]

 


 


 

 
ϵt ≤ ϵt ≤ ϵt ∀t ∈ [T ]

 


87
where: for all k ∈ N, d−k is equal to an historical value from the training sample, dt and d¯t
are lower and upper bounds on the demand, ϵt and ϵ̄t are lower and upper bounds on the
gaussian noise in period t, β0 is a constant such that β0 = µ− pk=1 βk µ, and the parameter
P

Γϵ ≥ 0 restricts the value of the noise vector ϵ. Notice that the uncertainty set restricts the
maximum cumulative deviation on the errors themselves: This ensures that each demand
sequence is consistent through time by linking the errors in the different time periods.

Remark 1. The interval [ϵt , ϵt ] can be any symmetric or asymmetric interval contain-
ing 0.

Remark 2. More advanced type of processes, such as the Autoregressive Integrated

Moving Average (ARIMA) or its seasonal version, involve to fit an ARMA process on
linear transformations of the original process. Embedding such structures thus only re-
quires simple modifications of the uncertainty set (5.13) without altering its nature (i.e.
bounded and convex polyhedron). For instance, if (Dt )t∈[T ] is modeled accurately by an
ARIM A(p, 1, q) process, it suffices to replace all terms of the form dj by dj − dj−1 in the
equality constraints.

3.1.2 Examples for special cases of demand processes

We focus on a simple case where the demand for a given product follows an AR(1)
process, or equivalently and ARM A(1, 0) process. The uncertainty set (5.13) based on
the AR(1) process uses the following demand relationship:

dt = β0 + β1 dt−1 + ϵt (5.14)

Assume the values of β0 , β1 and σϵ2 are fitted using their MLE on a training sample
containing N + 1 historical observations. Since ϵt has mean 0 for all t, we know from
equation (5.14) that
E(dt ) = β0 + β1 E(dt−1 ) (5.15)

and we can express the variance of the demand dt as:

1 − β12t
Var(dt ) = σϵ2 (5.16)
1 − β12

88
We use this formula to define the bounds dt and dt as follows:
v v
u 1 − β12t u 1 − β12t
u u
E(dt ) − Γd tσ 2
ϵ ≤ dt ≤ E(dt ) + Γd tσ 2
ϵ (5.17)
1 − β12 1 − β12

where Γd control the width of the interval. In what follows, we tune this parameter
empirically on the vector of demand [d−N , . . . , d0 ] in order to cover 95% of the possible
values of dt . Note that since parameter fitting assumes normally distributed errors, it is
quite natural to define a symmetric interval centered on the value E(dt ).
With ARMA(1,1), the induction becomes:

Dt = β0 + β1 Dt−1 + ϵt + θ1 ϵt−1 (5.18)

Hence, by induction, it is clear that one can express Dt as a linear combination of errors
ϵt and data d1 , . . . , dN . In particular, it is possible to derive a prediction interval for any
future period t > 0 with a given confidence level (1 − α)% using bootstrapping tech-
niques to obtain approximate confidence limits. This is well-documented in the literature
(see e.g. Brockwell and Davis (2016)) and readily implemented in numerous existing li-
braries dedicated to time series analysis (e.g. https://www.statsmodels.org/stable/
tsa.html in python or https://cran.r-project.org/package=tseries in R).

3.2 Joint demand and lead time uncertainty set

As mentioned in the introduction, we assume that the lead time may be influenced
by the value of past demands. Among others, Linear Regression (LR) is a simple way to
describe the dependency between multiple parameters. In our case, it consists in fitting
a linear model of the form: n
lt = ¯l +
X
bk dt−k + ζt (5.19)
k=0

where the uncertain lead time lt of period t is described by the constant value ¯l, the vector
of linear dependencies b and the random perturbation vector ζ lying in [ζ, ζ] and with
mean 0. Given a vector of N + 1 consecutive demand values [d−N , . . . , d0 ] and N + 1
consecutive lead time values [l−N , . . . , l0 ], one can use the Mean Squared Error (MSE)
estimator to fit the linear model and obtain the value of ¯l and b by minimizing the

89
function:
N 
1 X ˆlt − lt 2

M SE = (5.20)
N t=1

where lˆt is the prediction value defined by:

n
lˆt = ¯l +
X
bk dt−k (5.21)
k=0

Note that we can express the value of the lead time lt as :

T
X
lt = (1 − δtk ) (5.22)
k=t

We can then define a joint uncertainty set for demand and lead time by:
 


 δtk ≤ δtk+1 , ∀t = 1, . . . , T − 1, t ≤ k ≤ T (a) 



 

δtk ∈ {0, 1}, ∀t = 1, . . . , T, t ≤ k ≤ T (b) 

 

 


 

 


 T
X



Lmin ≤ (1 − δtk ) ≤ Lmax , ∀t = 1, . . . , T (c) 

 

 

 


 k=t 


 

 T n 

(1 − δtk ) = ¯l +
 X X 
bk dt−k + ζt , ∀t = 1, . . . , T (d) 

 

 

 
k=t k=0

 


 

UδLR = δ, d ζ ≤ ζt ≤ ζ, ∀t = 1, . . . , T (e)  (5.23)

 


 p
X q
X



dt = β0 + βk dt−k + θk ϵt−k + ϵt ∀t ∈ [T ] (f ) 

 

 

 
k=1 k=1

 


 

t √

 

 X 
| ϵj | ≤ Γϵ σϵ t ∀t = 1, . . . , T (g) 

 

 

 
j=1

 


 


 




 dt ≤ dt ≤ dt ∀t ∈ [T ] (h) 



 

 

ϵt ≤ ϵt ≤ ϵt ∀t ∈ [T ] (i) 

where constraints (a) ensure that the quantity received until period t + 1 is greater
than or equal to the quantity received until period t. Constraints (b) forces the integrality
of δ, which in turn implies the integrality of lead times. Constraints (c) and (d) bound
the value of lead times and constraints (e) bound the possible values of the perturbations
ζt . Note that this set of constraint allows for order crossover, (i.e. we can receive the
order qt before precedent orders qt−n for all n ∈ [1, Lmax − 1]). Finally, constraints (f)-(i)
are constraints related to the demand.

90
3.3 A constraint generation approach

In uncertainty set (5.12), the binary constraint on δ is relaxed to define a convex

uncertainty set. This property is required in classic robust optimization to use the dual
formulation of the subproblem. In the uncertainty set (5.23), the constraint (d) integrates
the correlation between demand and lead time. In order to use the linear model (5.19)
obtained by minimizing function (5.20) we need to ensure the integrality of the lead time
values (constraint (b)). As a result, the proposed data-driven uncertainty set is a finite
discrete set and solving the problem using a dual reformulation approach is not possible.
However, one can adapt the constraint generation approach introduced in Bienstock and
Ozbay (2008) and used in Thorsen and Yao (2017) to solve the robust formulation of
the problem with the joint discrete uncertainty set. This approach rely on 1. a master
problem that take decisions according to a set of uncertain scenarios and 2. a sub problem
generating worst case scenarios to enrich the master problem. Both problems are solved
iteratively by enriching a set Ω of lead time and demand scenarios until an optimal robust
solution is found. For a given set of scenarios Ω, let

T ×|Ω|
n o
X (Ω) = (x, q, y) : x ∈ {0, 1}T , q ∈ RT+ , y ∈ R+

, the master problem of the inventory planning model under uncertain lead time and
demand is given by:

min z (5.24)
x,q,y,z
X
s.t. z ≥ sxt + cqt + ytω ∀ω ∈ Ω (5.25)
t∈[T ]
t
!
(δitω qi dωi )
X
ytω ≥ h I0 + − ∀t ∈ [T ] , ∀ω ∈ Ω (5.26)
i=1
t
!
(δitω qi dωi )
X
ytω ≥ −b I0 + − ∀t ∈ [T ] , ∀ω ∈ Ω (5.27)
i=1

qt ≤ M x t ∀t ∈ [T ] (5.28)
(x, q, y) ∈ X (Ω) (5.29)

By solving this problem, we obtain the optimal decision vectors x∗ and q ∗ that are used
to identify a new worst-case scenario to add to Ω. This scenario is obtained as the optimal

91
solution to the following MILP:

T
sx∗t + cqt∗ + Ht + Bt
X
max (5.30)
δ,d,H,B,p
t=1
t
!
(δit qi∗
X
s.t. Ht ≥ h I0 + − di ) ∀t ∈ [T ] (5.31)
i=1
t
!
(δit qi∗
X
Ht ≤ h I0 + − di ) + Mt (1 − pt ) ∀t ∈ [T ] (5.32)
i=1

Ht ≤ Mt pt ∀t ∈ [T ] (5.33)
t
!
(δit qi∗ − di )
X
Bt ≥ −b I0 + ∀t ∈ [T ] (5.34)
i=1
t
!
(δit qi∗
X
Bt ≤ −b I0 + − di ) + Mt pt ∀t ∈ [T ] (5.35)
i=1

Bt ≤ Mt (1 − pt ) ∀t ∈ [T ] (5.36)
Ht , Bt ≥ 0 ∀t ∈ [T ] (5.37)
pt ∈ {0, 1} ∀t ∈ [T ] (5.38)
δ, d ∈ UδLR (5.39)

where Mt are sufficiently large values. This sub-problem maximize the inventory costs
obtained with the optimal decisions x∗ , q ∗ of the master problem by selecting the worst
realization of demand and lead time in the uncertainty set.
The procedure for this constraint generation approach is given by Algorithm 1. Bi-
enstock and Ozbay (2008) proved that this algorithm converges to an optimal robust
solution in a finite number of iterations when the uncertainty sets are non-empty and
bounded.

4 Computational experiments
The objective of this section is to evaluate the solutions provided by the robust models
introduced up to now and compare their performances against a sample average approach
(presented in Section 4.1). In order to evaluate the impact of each contribution, we
consider two different scenarios set of experiments. In the first set of experiments -
Section 4.2), we consider that the demand follows an autoregressive process AR(1) while
lead times are uniformly distributed in an interval. In the second set of experiments

92
Algorithm 1 Adversarial Approach Bienstock and Ozbay (2008)
1. Initialization:
Select an initial set of scenarios Ω
Define tolerance value ϵ
Set LB = 0 and U B = ∞
2. Master problem:
Solve the master problem with scenario set Ω to obtain x∗ , q ∗ and Z ∗
Set LB = Z ∗
3. Sub problem:
Solve adversarial sub-problem with x∗ and q ∗ to obtain δ ∗ , d∗
Set U B = min(U B, Objective value of sub-problem)
4. If U B − LB > ϵ:
Add scenario (δ ∗ ,d∗ ) to Ω, and return to Step 2.
5. If U B − LB ≤ ϵ:
Terminate and return q ∗ .

(Section 4.3), we study the impact of potential correlations between demand and lead
time.

4.1 Sample average approach

The Sample Average Approximation (SAA) model is formulated as a two-stage model
in which the first stage decisions defines the ordering periods x as well as the quantities
ordered q, while the second stage decisions computes the holding and backlogging costs
y for each scenario. The model is given by the following MIP:

T
1 XX
min sxt + cqt + ytω (5.40)
x,q,y |Ω| ω∈Ω t=1
t
!
(δitω qi dωi )
X
s.t. ytω ≥ h I0 + − ∀t ∈ [T ] , ∀ω ∈ Ω (5.41)
i=1
t
!
(δitω qi dωi )
X
ytω ≥ −b I0 + − ∀t ∈ [T ] , ∀ω ∈ Ω (5.42)
i=1

qt ≤ M x t ∀t ∈ [T ] (5.43)
(x, q, y) ∈ X (Ω) (5.44)

in which Ω is the set of demand and lead time scenarios. The objective function (5.40)
aims at minimizing the average cost over the set of scenarios. Constraints (5.41) and
(5.42) compute the costs due to storage and backlogging in each period and each scenario.
Finally, constraints (5.43) link the quantity ordered qt and the binary variable xt .
Clearly, the quality of an SAA solution depends heavily on the set of scenarios Ω that

93
Table 5.1 – Costs obtained by SAA solutions with AR(1) and Normal demand scenarios
for instances of 20 periods without fixed cost, ratio hb = 8, σϵ = 10 when β1 vary.

Average Median Standard deviation Q0.99

β1
Normal AR(1) Normal AR(1) Normal AR(1) Normal AR(1)
0.3 695 639 551 573 488 308 2892 2055
0.5 1092 780 702 700 982 414 4955 2624
0.8 2916 1325 1631 1154 2843 864 12219 5198

is included in the formulation. In particular when the demand follows a specific process
such as the ARMA described before, one may improve the performances of the formulation
by generating trajectories using the information embedded in the historical data . As a
matter of example, Table 5.1 shows the costs of solutions obtained with SAA models
on an instance with 20 periods, no fixed cost and a ratio hb = 8 in which the demand
follows an AR(1) process. We study three types of instances corresponding to different
degree of auto-correlation and each line present the average values obtained over a set of
10 instances. For each type, we generate a scenario set containing 500 scenarios whose
demand distribution follows either a Normal distribution or an AR(1). Each parameters
used to generate scenarios (e.g. β1 , mean value, standard deviation) are estimated on
historical data sets of 100 demand realizations. The solutions found by each formulation
are then evaluated on 10 000 different trajectories and for each value of β1 the process is
repeated on 50 instances.
The results indicate that the SAA using scenarios generated with the AR(1) process
lead to better solutions than the one based on an iid demand generation, both on average
and in the worst case. For example when β1 = 0.8 the average and extreme (0.99 quantile)
costs obtained with the normally distributed training sample are more than twice the ones
obtained with the one generated based on an autoregressive process.

4.2 Uncertain correlated demand

The results from Table 5.1 suggest that being able to precisely characterize the demand
process leads to a better scenario generation. However in practice, the only information
available to define an appropriate model and tune its parameters consists of the historical
demand values. Even when the nature of the undelying process is known, biased estimates
of the parameters may affect scenario generation and deteriorate the performances of the

94
 Table 5.2 – Instance configurations for the first set of experiments

Name Definition Value

T time horizon 20
c unit ordering cost 0.1
h unit holding cost 0.3
b unit backlogging cost 2.4
s fixed ordering costs {0, 45, 105, 165}
T BO time between orders {0, 1, 2, 3}
β1 autocorrelation factor {0.3, 0.5, 0.8}
β0 constant factor 100(1 − β1 ) (gives an average demand value of 100)
σϵ standard deviation of the perturbation 10

SAA. In this context, robust optimization may be seen as an alternative to prevent decision
makers from paying extreme inventory costs because of these errors. In the remainder of
this section, we assume that all the models use a common data set containing N = 100
consecutive demand (and leadtimes in Section 4.3) observations. In particular when the
models have access to the specific demand process, they define its parameters as their
MLE using the data available (and their possible bias). On each instance, we compare
the performances of the SAA formulation against two robust models that use either a
budget-based or an ARMA-based uncertainty set.
In this first set of experiment, the demand evolves according to an AR(1) process with
several combination of parameters as described in Table 5.2. Note that the different values
of fixed cost s are computed to control the Time Between Order (TBO), which corresponds
to the expected number of periods between two consecutive orders in order to minimize the
total costs. For each combination, we generated 50 distinct instances, i.e. LSP that share
the same parameters except their training sample of historical demands. For each instance,
we compare the solution obtained by the SAA described in Section 4.1 a Sample Average
Approximation (SAA), a robust approach with the classic uncertainty set UdΓ described
ARM A(1,0) AR(1)
in (5.11) and the AR(1) based uncertainty set described above Ud = Ud .
We evaluate the solutions obtained by each method on an independent test set of
K = 105 scenarios generated using Monte Carlo simulation and we report the average,
0.95- and 0.99-quantiles of the empirical cost achieved by each solution over the test
sample. More precisely, we report the absolute cost evaluated by simulation for the SAA
model while for both robust approaches we report the relative cost difference from the
SAA solution in %. For the classical robust model, the budget parameters Γt are set for all
t ∈ [T ] with the objective to include the cumulative demand in the uncertainty set with

95
 AR(1)
Table 5.3 – Costs obtained with SAA, UdΓ and Ud for instances with fixed cost s = 0
when β1 vary.

Average Q0.95 Q0.99

β1
AR(1) AR(1) AR(1)
SAA UdΓ Ud SAA UdΓ Ud SAA UdΓ Ud
0.3 690.28 +33.44% +0.63% 1385.51 -15.46% -16.48% 2266.19 -29.39% -23.97%
0.5 859.47 +20.16% +3.57% 1819.01 -15.99% -16.3% 2973.26 -23.34% -21.3%
0.8 1480.14 +11.83% +10.29% 3473.42 -6.28% -20.26% 5875.9 -10.64% -27.23%

an approximate probability of 0.95, under the assumption of IID normally distributed

demands. Assuming that the training sample empirical mean and standard deviation
are d¯ and σ̂d , respectively, we denote for all t ∈ [T ] and k = t, . . . , N + 1: ∆k (t) =
Pt ¯ ¯
j=1 |d−k+j − d|, the cumulative absolute deviation from d of a sequence of t consecutive
demands starting from index −k in the training data set. We then define for all t ∈ [T ]
the budget value Γt = min{γ ≥ 0 : N k=t 1{∆k (t)≤2γ σ̂d } ≥ 0.95(N + 1 − t)}. For the AR(1)-
P

based uncertainty set, we use the normality assumption on the error terms and define
P √
Γϵ = 2 so that P | tk=1 ϵt | ≤ Γϵ σϵ t ≈ 0.95.
Table 5.3 shows the costs of solutions obtained with different degree of autocorrelation
β1 when there is no setup costs. In this situation, we observe that the SAA model always
provide the lowest average cost but also the highest extreme costs. For example, when
β1 = 0.3, the 0.99-quantile of the SAA solution is more than three times its average
cost. In comparison, the robust approaches provides much higher average costs but with
a significant reduction of extreme costs. Using the AR(1)-based uncertainty set provides
the lowest average costs and 0.95-quantile among the two robust approaches and across
all values of β1 . Interestingly, the classic robust approach outperforms the AR(1)-based
in case of weakly auto-correlated processes (β1 ∈ {0.3, 0.5}) but this phenomenon fades
as β1 increases. In particular when β1 = 0.8 the solutions derived using the AR(1)-based
approach give the best 0.99-quantile performance. Generally speaking, it appears that the
budget-based uncertainty set UdΓ cannot lead to both robust solutions and good average
AR(1)
performances simultaneously compared to Ud .
Table 5.4 summarizes the costs obtained with different values of setup cost s when
the autocorrelation factor β1 = 0.8. We observe that increasing setup costs decrease
the cost difference across all methods without modifying the previous observations. The
SAA model still provide the lower average cost and still lead to very high 0.95- and 0.99-
quantile. For example, when the setup cost s = 45, the AR(1) based robust solution

96
 AR(1)
Table 5.4 – Costs obtained with SAA, UdΓ and Ud for instances with β1 = 0.8 when
the fixed cost s vary.

Average Q0.95 Q0.99

s
AR(1) AR(1) AR(1)
SAA UdΓ Ud SAA UdΓ Ud SAA UdΓ Ud
0 1480.14 11.83% 10.29% 3473.42 -6.28% -20.26% 5875.9 -10.64% -27.23%
45 2022.89 4.54% 2.17% 3665.55 -9.14% -14.58% 5785.06 -16.5% -22.87%
105 2517.81 2.57% -0.09% 4151.87 -5.7% -11.95% 5980.36 -9.99% -18.53%
165 2825.64 2.43% 0.6% 4304.56 -6.1% -11.32% 5970.0 -11.15% -19.74%

increase the average cost by 2.17% while it reduces the average 0.99-quantile by 22.87%.

In order to compare in more details the relative performances of the three methods, we
compare the empirical superquantile of the cost achieved by each method over the test set.
Superquantile is a well-established measure to evaluate the risk taken by a decision maker.
Also known as conditional value at risk or expected shortfall in portfolio optimization, it
corresponds to the average cost above a given quantile qα of the cost distribution. For a
given random variable X, the superquantile is given by:

1
 h i 
Sα (X) = E [X|X ≥ qα (U )] = min η + E (X − η)+
η∈R 1−α

, with its empirical counterpart for a sample x̂ = (x̂1 , . . . , x̂K ) of size K:

K
( )
1
(x̂i − η)+
X
Ŝα (x̂) = min η +
η∈R (1 − α)K i=1

In the later, we compare the costs obtained by the SAA approach and the classic bud-
get based robust approach with the AR(1)-based robust approach by computing the
relative gap between the values of their respective superquantile for various values α ∈
[0, 1]. Recall that for each combination of parameters, we create 50 historical trajectories
(i) (i)
[d−N , . . . , d0 ] that are used to obtain one solution to the corresponding LSP for each
method. Each solution computed for instance i is then evaluated on K distinct out-of-
sample trajectories of the process, giving three vectors containing their own realized costs
(i) (i,1) (i,K) (i) (i,1) (i,K) (i) (i,1) (i,K)
ẑSAA = (ẑSAA , . . . , ẑSAA ), ẑΓ = (ẑΓ , . . . , ẑΓ ) and ẑAR(1) = (ẑAR(1) , . . . , ẑAR(1) ). For a
given instance i = 1, . . . , 50, one can thus compare the relative superquantile of the three

97
methods by computing
   
(i) (i)
(i)
Ŝα ẑSAA − Ŝα ẑAR(1)
ρSAA (α) = 
(i)
 (5.45)
Ŝα ẑAR(1)
   
(i) (i)
(i)
Ŝα ẑΓ − Ŝα ẑAR(1)
ρΓ (α) = 
(i)
 (5.46)
Ŝα ẑAR(1)

Notice that we take as a reference the cost returned by the AR(1)-based robust formulation
to facilitate the evaluation of its performances compared to the other two. In Figure 5.1
and 5.2, we represent the box plot across 50 instances of the relative gaps defined in (5.45)
and (5.46) on the superquantiles of the different methods when α increase. Figure 5.1
shows the evolution of this ratio for different values of β1 with no setup costs and all
other parameters fixed, while Figure 5.2 compares the same evolution for a fixed β1 as
the setup cost increases. We observe in line with Table 5.3 and 5.4 that SAA provides
the best average cost in most instances while the classic robust approach provide good
performances on extreme values of the cost distribution (α = 0.9). Another important
observation is that as β1 increases, the cost difference between all methods decreases (on
average over all instances) as stated by table 5.3. However both the SAA and the budget-
based robust formulation compute solutions with much stronger positive deviations from
the cost achieved by the AR(1)-based robust solution, while improvements over the latter
by one of the former is both rare and of less magnitude across all values of α and all the
combinations tested.

4.3 Uncertain demand-dependent lead times

In this second set of experiments, we consider uncertain demand and lead times and
assume that a positive correlation exists between these two types parameters. In the
remainder, we denote by UδΓ the use of classical uncertainty sets presented in Section 2 for
demand and lead time while UδLR refer to the data-driven joint uncertainty set introduced
in Section 3.

Lead time generation Given a set of N + 1 consecutive demand [d−N , . . . , d0 ] we

generate a set of N + 1 consecutive lead time values [l−N , . . . , l0 ] taking values in [L] with
the following formula:
lt = min{L, max(1, Zt2 ∗ L}} (5.47)

98
 (a) s = 0, β1 = 0.3

(b) s = 0, β1 = 0.5

(c) s = 0, β1 = 0.8

Figure 5.1 – Evolution of the relative Ŝα as α increase for solutions obtained with SAA
and UdΓ for an instance with ratio hb = 8, σϵ = 10 and s = 0 for different values of β

99
 (a) s = 45, β1 = 0.8

(b) s = 105, β1 = 0.8

(c) s = 165, β1 = 0.8

Figure 5.2 – Evolution of the relative Ŝα as α increase for solutions obtained with SAA
and UdΓ for an instance with ratio hb = 8, σϵ = 10 and β1 = 0.8 for different values of setup
cost s

100
Table 5.5 – Costs obtained with SAA, UδΓ and UdLR for instances with ratio hb = 8, σϵ = 10,
β1 = 0.5 when the fixed cost s = 105 and the maximum lead time L varies.

Average Q0.95 Q0.99

L
SAA UδΓ UδLR SAA UδΓ UδLR SAA UdΓ UδLR
1 1952.36 +7.71% +1.41% 2500.95 -10.79% -12.54% 3044.97 -18.04% -19.02%
3 2935.96 +6.74% +3.02% 4097.26 -7.61% -10.04% 5351.84 -15.77% -16.86%
5 3002.59 +4.87% +3.13% 4132.72 -7.28% -9.8% 5710.53 -13.79% -13.85%

( nk=1 dt−k )−D

P
, where Zt = D−D
and D and D are respectively the maximum and minimum
value of the sum of n consecutive demands nk=1 dt−k , evaluated on the training sample
P

of N + 1 demands. For the current set of experiment, we set n = 3 periods which mean
that the lead time of order depends on the level of the demand over the three previous
periods.

Numerical results Table 5.5 shows the average cost obtained by the three methods
applied on 50 instances sharing the same configuration. We consider a set of instances
with β1 = 0.5, s = 105 and for which the maximum lead time L vary from 1 to 5. The
results in this table remain consistent with the previous observations as SAA formulations
lead to the lowest average cost and the highest extreme costs. The average cost difference
between the SAA and the classic robust solutions decreases when the maximum lead time
increases, while the difference between SAA and the solutions provided by our method
increases.
The joint demand and lead time uncertainty set presents clear advantages over the
one used by budget-based robust model, as it is more robust and performs better on
average than the latter. However in our set of experiments, the cost difference between
both robust models decreases as the maximum lead time L increases. In particular when
L = 5 the average cost difference is less than 2% between both robust solutions and the
0.95-quantiles of costs are separated by less than 3%.

5 Conclusion and perspectives

In this chapter, we study an inventory management problem in which the demand level
and the lead time of orders are uncertain and can be correlated. We introduce a paramet-
ric approach to build uncertainty set in which demand is assumed to follow a stochastic

101
process. We use classical autoregressive and moving average models to formulate uncer-
tainty sets that integrates information on the demand process such as autocorrelation.
We propose numerical experiments in which the demand follows an AR(1) process and
show that our approach exploiting relationships between consecutive demands outper-
forms a classic robust approach inspired Bertsimas and Sim (2003, 2004) applied directly
to the demand values. When compared to a standard sample average approximation
scheme, it significantly reduces the risk to incur an extremely high cost at the expense of
a mild increase the average cost of solutions.
Our results suggest that lead time uncertainties affect negatively the competitiveness
of using our AR(1) based uncertainty set for demand uncertainty, in particular when the
maximum lead time is high. Integrating a potential correlation between demand and
lead time does not appear to provide a particular advantage in comparison to a classic
budget-based uncertainty set when positive correlation is considered.
With this first study, we show that using statistical model such as ARMA (or, by
extension, ARIMA) to integrate correlation information improves the performance of
robust models when uncertain parameters are correlated. Several interesting extension of
our work should be considered in future research. First of all, we restricted our experiment
to an AR(1) model but our framework can be applied to more general ARIMA models
which may further advantage structured uncertainty sets over generic ones. An other
direct extension of our work would be to consider the application of such uncertainty sets
in a dynamic context in which ordering quantity are adapted as the planning horizon unroll
and some uncertain parameters are realized. Finally, the application of this approach on
a real case study would be particularly interesting to conclude on the practical relevance
of the method in an industrial context.

102
 Chapter 6

C ONCLUSIONS AND P ERSPECTIVES

In this thesis we study the application of data-driven robust models on several opti-
mization problems arising in supply chain management. The literature review presented
in Chapter 1 highlights that this paradigm is still young and despite promising results in
various application domains, the use of such data-driven models in supply chain manage-
ment and production planning remains relatively scarce. The main difficulty that prevents
these new data-driven robust models to translate into methods that are applicable in prac-
tice comes from their lack of scalability. Indeed data-driven robust optimization may not
be a viable option in SCM considering the (potentially large) MIP formulations of the
optimization problems that arise in this context.
In this manuscript, we demonstrate that using data-driven uncertainty sets to bound
the behaviour of uncertain parameters of optimization problems more precisely can pro-
vide significant reduction of the conservatism of robust solutions to SCM problems. In
addition, we address the scalability issue through the introduction of approximate data-
driven approaches, with the objective to improve the tractability of the resulting robust
optimization models. This section recalls the main conclusions drawn from our work,
both in terms of scientific contributions and future research perspectives susceptible to
turn these techniques into useful tools for industrial applications.

Chapter 3. Our study of the SVC-based uncertainty set exhibits the lack of scalability
of this approach, despite its interest for various applications (see Chapter 1, Section 1.3).
We identify the large number of addition variables and constraints introduced in SVC-
based robust formulations as one of the main causes of this limitation and propose an
alternative approximate model that reduces the size of the MILP formulations while
maintaining a good-quality representation of uncertain parameters. We compare the
solutions obtained by the SVC-based approach and our approximation on three distinct
optimization problem (MILP) in order to demonstrate the efficiency of our approach.
Our computational results show that our approximation succeeds in producing solu-

103
tions that are comparable to the solution obtained with the SVC-based uncertainty set.
For example, on a robust knapsack problem with uncertain correlated weights, our ap-
proximation provide solutions which objective values are less than 1% higher than the
solutions of the SVC-based uncertainty set on most of the instances. On the other hand,
the computation time needed to solve these new robust models is reduced by more than
90%. In addition, our approximation scheme uses a single parameter that scales with
the dimension of uncertain parameters and efficiently controls the trade-off between the
quality of the approximation and the computation time (i.e. the size of the robust coun-
terpart).
As a result, this study extends the range of application of the SVC-based approach
to more complex and larger optimization problems. In addition, as it rely on the same
formulation as the initial uncertainty set, it is likely to be compatible with any other
improvement of the SVC-based uncertainty set that have been proposed in the literature
(see Chapter 1 Section 1.3 for examples).

Chapter 4 This chapter provides a practical application of the approximate SVC-based

uncertainty set. We focus on a problem faced by a manufacturer who must plan the
assembly of multiple products in order to meet a known demand over a finite planning
horizon. The picking and delivery of the components necessary to obtain the end items
are subcontracted to an external logistic service provider. The goal is to define the orders
in order to respect the maximum time capacity dedicated to the corresponding picking
operations when their duration is uncertain.
We propose a MIP formulation and develop two robust models that integrate un-
certainties to solve this problem. The first one is based on the budget-based approach
of Bertsimas and Sim (2003, 2004) and allows to reduce the impact of uncertain pick-
ing times on the assembly line by providing solutions with a higher demand satisfaction
compared to a deterministic solution. As an alternative robust solution approach, we
apply the approximate SVC-based one introduced in Chapter 3. A set of experiments
on instances with up to 10 products and 20 components shows that our approximate
SVC-based uncertainty set outperforms the classic budget based approach on this prob-
lem. Specifically, it simultaneously reduces significantly the risk of backlogging demand
and leads to a lower number of components held on the border of the assembly line. For
example, on an instance with 5 components and 10 products over a planning horizon of
14 periods, our approach can reduce the average and extreme (0.95-quantile) total cost

104
incurred by producers by 20%.

Chapter 5 In the first two contributions, we make use and extend a data-driven robust
model based on a non-parametric approach (i.e. no assumptions are made on the uncertain
parameters). However, this type of approach requires a large amount of historical data to
be efficient, which restricts their range of application in supply chain management. For
example when the data consists of time-dependent uncertainties, e.g. in the form of a
demand process, gathering large and reliable data set is often impossible as the process
that drives the parameters may evolve over time (e.g. new products on the market, crisis,
...). In this context, we introduce a novel uncertainty set built upon the assumption that
uncertain parameters are driven by a specific (but unknown) stochastic process that is
commonly used in predictive approach.
To illustrate the efficiency of our approach, we study a simple single-item lot sizing
problem under uncertain and autocorrelated demand and an extension of this problem in-
tegrating an uncertain, demand-dependent lead-time. In addition, we rely on a dedicated
cutting plane method that is known to be less conservative than the dual reformulation
approach to solve this problem (Thorsen and Yao (2017)).
Our experiments show that our uncertainty set results in a significant reduction of
conservatism compared to a classic budget-based approach. The solutions that are found
using our parametric data-driven approach are of comparable or lower extreme costs with
a significant reduction of cost on average. We also compare the solutions obtained to
those of a sample average approximation method whose scenarios integrate the same
correlation information that our robust model. The experimental results show that the
sample average approach provides the best performance on average but that it performs
poorly in extreme cases. In addition, when the correlation is high, the performances of
the sample average approach display more variability than our robust formulation. For
some instances, the solution of the latter is better than the former both on average and
in extreme cases.
In a second set of experiments, we integrate uncertain and demand-dependent lead
time to the problem. The results show that this type of uncertainty significantly impacts
the costs of the solutions obtained with all three methods, which all adopt a « robust »
policy with regards to uncertain lead times. The conclusions drawn from the first set
of experiments remain unchanged when we couple uncertain lead times with uncertain
correlated demand.

105
Perspectives Scalability of DDRO approaches: The scalability issue is a major draw-
back of most data-driven robust optimization methods. We show in this thesis that
approximating uncertainty sets may address this problem and extend the range of appli-
cation of existing methods. While our proposition is only applicable to a specific type of
uncertainty sets, we believe that the core idea of our method may be adapted to other
families of uncertainty sets. Extending this framework may improve the range of appli-
cations for which robust counterparts remain tractable in practical settings. In addition,
most existing studies on similar data-driven methods rely on the dual reformulation to
obtain solutions of robust models and the existing works comparing dual reformulation
and cutting plane approaches to solve robust models is restricted to classic uncertainty
sets (i.e. budget based, ellipsoidal). Moreover, their is no consensus on which method
should be used and the efficiency of both methods seems to vary among uncertainty sets
and optimization problems (see Chapter 1 Section 1.1). Thus, solution methods based on
cutting plane approaches could help in order to improve the scalability of DDRO methods.
As a small step on this path, we studied the difference of computation time to solve the
robust binary knapsack problem of Chapter 1 with dual reformulation and with a cutting
plane approach as described in Chapter 1 Section 1.1. Figure 6.1 shows the computation
time needed to solve the binary knapsack problem with both methods. We present results
obtained with the SVC-based uncertainty set (6.1a) and our approximation (6.1b). In
Figure 6.1a each boxplot correspond to a set of 10 instances and the computation time
evolve with the parameter ν of the SVC based approach. In Figure 6.1b each boxplot
correspond to a set of 90 problems instances and the computation time evolves with the
parameter K of our approximation approach. This first set of experimental results clearly
shows that in this context, coupling the SVC-based uncertainty set or our approximation
with a cutting plane approach greatly reduces the computation time needed to solve the
resulting formulations.
Improved cutting plane algorithm for LSP: As described in Chapter 1 and Chapter 5,
the dedicated cutting plane algorithm introduced by Bienstock and Ozbay (2008) for
lot sizing problems solves a relaxed robust problem with only a subset of realizations
selected from their uncertainty sets. At each iteration, a new scenario is added to the
formulation in order to converge to the optimal solution of the true min-max robust
problem. One limitation of this procedure is that given a subset of scenarios, there could
exist multiple optimal solutions to the min-max problem, each with a different average
objective value over the subset of scenarios. As an alternative to this pure « worst case »

106
 (a) SVC-based uncertainty set (b) Approximated uncertainty set

Figure 6.1 – Computation time of dual reformulation and cutting plane method for the
binary knapsack problem with 30 items using the SVC-based uncertainty set and our
approximation.

perspective, one could think of a bi-objective approach that focuses on the worst case
performance as its first objective, then on the average cost over the subset of scenarios
selected in a lexicographic order. The obtained solution would then have the same worst
case performance but likely be less conservative on average. To go even further, one could
also think of implementing this method into more general two stages robust problems
solved with a column and constraint generation approach such as the one presented in
Zeng (2011).
Adjustable data-driven robust optimization: This thesis focuss on static decision con-
text in which all decisions are fixed prior to the realization of uncertain parameters.
During the last two decades, adjustable robust models have been proposed to obtain so-
lutions that depend on the realization of uncertain parameters (see Chapter 1 Section
2.3). Such solutions may be of interest for many applications in SCM and are known to
be less conservative than static decisions making. Studies coupling such approaches with
data-driven methods are still scarce and often applied to other fields. A direct extension
of our work would then be to consider adjustable versions of the problem introduced in
this thesis, e.g. to tackle the LSP from Chapter 5. For instance, a direct extension that
consists in fixing the setup before the beginning of the planning horizon and adjusting
the lot size decisions dynamically as the uncertainty is revealed would likely improve the
solution quality drastically. In addition for this specific case, we expect that the more
precise representation of the uncertainty developed in Chapter 5 would favor the dynamic
counterparts of these static robust lot-sizing models.

107
 Link with constraint satisfaction level: Most data-driven robust methods aims at defin-
ing an uncertainty set U(D) that covers a given portion of historical data D. While this
may give an accurate estimation of the probability that the realization of random param-
eters falls within such sets, this portion is often an underestimate of the probability of
satisfaction for the corresponding robust constraint. A more interesting property would be
to design sets that depend on both the uncertain parameters and the decision variables
in order to define a probabilistic guarantee of constraint satisfaction instead. In other
words, considering that the set D describe the true distribution of uncertain parameters,
a desirable property of a DDRO approach would be of the form:

P(a ∈ U(D, x)) ≥ P(a⊺ x ≤ b) (6.1)

Existing works rely on the concept of robustness in distribution to derive ambiguity

sets that allow to bound this probability, but they usually result computationally heavy
models. In practice it is still difficult to estimate the true constraint satisfaction level
obtained from a particular uncertainty set. Whether it is possible to design a data-driven
uncertainty set that provides a tight probabilistic bound on a given constraint satisfaction
thus remains an open research question. Such bounds would be of particular interest for
practitioners in order to tune data-driven uncertainty set to their specific needs and avoid
overconservative decisions.

108
 A PPENDICES

A Support Vector Clustering and Kernel function

In this section, we detail the formulation of the Quadratic Problem solved in the
Support Vector Clustering algorithm and the kernel function introduced by Shang et al.
(2017). The SVC algorithm seeks to minimize the radius R of a sphere enclosing a portion
1−ν of the available data {ϕ(u(i) )} where ϕ(.) is a mapping to a feature space F in higher
dimension. The corresponding optimization problem is given by:

N
1 X
min R2 + ξi
N ν i=1
s.t. ||ϕ(u(i) ) − p||2 ≤ R2 + ξi u(i) ∈ D
R∈R
p ∈ Rm
ξi ≥ 0 i = 1, . . . , N

where p is the center of the sphere and ξ is a vector of slack variables defining if a
given data-point u(i) ∈ D resides into the sphere (i.e. if ξi > 0 then u(i) is outside of
the sphere). Introducing lagrangian multipliers α lead to the following equivalent dual
quadratic problem:

N X
N N
αi αj K(u(i) , u(j) ) − αi K(u(i) , u(i) )
X X
min
i=1 j=1 i=1

s.t. 0 ≤ αi ≤ 1/N ν i = 1, . . . , N
N
X
αi = 1
i=1

where K(u(i) , u(j) ) = ϕ(u(i) )T ϕ(u(i) ) is the kernel function. In order to obtain a poly-
hedral description of the cluster, Shang et al. (2017) introduced the Weighted General

109
Intersection Kernel given by:
m
(i) (j)
lk − ∥Q(u(i) − u(j) )∥1
X
K(u , u ) = (2)
k=1

1
where Q is called the weighting matrix and is defined by Q = Σ− 2 with Σ the covariance
matrix and where lk is such that it satisfy the constraint lk > max qkT u(i) − min qkT u(i) .
1≤i≤N 1≤i≤N
From the optimal solution α∗ of this QP we identify a set of support vectors S =
n o n o
u(i) |αi∗ > 0 and a set of boundary support vectors B = u(i) | 0 < αi∗ < 1/N ν that
are used to define the SVC-based uncertainty set as follows:
( )
αi∗ ∥Q(u − u(i) )∥1 ≤ θ
X
USVC = u
i∈S

αi∗ ∥Q(u(j) − u(i) )∥1 .

P
, where θ = maxj∈B i∈S

110
 B IBLIOGRAPHY

Abdel-Aal, M. (2019). A robust capacitated lot sizing problem with setup times
and overtime decisions with backordering allowed under demand uncertainty. IFAC-
PapersOnLine, 52:589–594.

Agra, A., Poss, M., and Santos, M. (2018). Optimizing make-to-stock policies through a
robust lot-sizing model. International Journal of Production Economics, 200:302–310.

Aharon, B.-T., Boaz, G., and Shimrit, S. (2009). Robust multi-echelon multi-period
inventory control. European Journal of Operational Research, 199:922–935.

Alem, D., Curcio, E., Amorim, P., and Almada-Lobo, B. (2018). A computational study
of the general lot-sizing and scheduling model under demand uncertainty via robust
and stochastic approaches. Computers & Operations Research, 90:125–141.

Aouam, T. and Brahimi, N. (2013). Integrated production planning and order acceptance
under uncertainty: A robust optimization approach. European Journal of Operational
Research, 228(3):504–515.

Atan, Z., Ahmadi, T., Stegehuis, C., de Kok, T., and Adan, I. (2017). Assemble-to-order
systems: A review. European Journal of Operational Research, 261(3):866–879.

Bandi, C. and Bertsimas, D. (2012). Tractable stochastic analysis in high dimensions via
robust optimization. Math. Program, 134:23–70.

Barlas, Y. and Gunduz, B. (2011). Demand forecasting and sharing strategies to reduce
fluctuations and the bullwhip effect in supply chains. The Journal of the Operational
Research Society, 62:458–473.

Ben Ammar, O., Dolgui, A., Hnaien, F., and Louly, M. A. (2013). Supply planning and
inventory control under lead time uncertainty: A review. IFAC Proceedings Volumes,
46(9):359–370.

111
Ben-Tal, A., Golany, B., Nemirovski, A., and Vial, J.-P. (2005). Retailer-supplier flexible
commitments contracts: A robust optimization approach. Manufacturing & Service
Operations Management, 7:248–271.

Ben-Tal, A., Goryashko, A., and Guslitzer, E. Nemirovski, A. (2004). Adjustable robust
solutions of uncertain linear programs. Mathematical Programming, 99:351–376.

Ben-Tal, A., Hazan, E., Koren, T., and Mannor, S. (2015). Oracle-based robust optimiza-
tion via online learning. Operations Research, 63:628–638.

Ben-Tal, A. and Nemirovski, A. (1998). Robust convex optimization. Mathematics of

Operations Research, 23(4):769–805.

Ben-Tal, A. and Nemirovski, A. (2000). Robust solutions of linear programming problems

contaminated with uncertain data. Math. Programming, 88:411–424.

Bertsimas, D., Dunning, I., and Lubin, M. (2016). Reformulation versus cutting-planes
for robust optimization. Computational Management Science, 13:195–217.

Bertsimas, D. and Georghiou, A. (2015). Design of near optimal decision rules in multi-
stage adaptive mixed-integer optimization. Operations Research, 63(3):610–627.

Bertsimas, D., Gupta, V., and Kallus, N. (2017). Data-driven robust optimization. Math-
ematical Programming, 167:235–292.

Bertsimas, D., Pachamanova, D., and Sim, M. (2004). Robust linear optimization under
general norms. Operations Research Letters, 32:510–516.

Bertsimas, D. and Sim, M. (2003). Robust discrete optimization and network flows.
Mathematical programming, 98(1-3):49–71.

Bertsimas, D. and Sim, M. (2004). The price of robustness. Operations Research, 52:35–53.

Bertsimas, D. and Thiele, A. (2004). A robust optimization approach to supply chain

management. In Integer Programming and Combinatorial Optimization, pages 86–100,
Berlin, Heidelberg. Springer Berlin Heidelberg.

Bertsimas, D. and Thiele, A. (2006). A robust optimization approach to inventory theory.

Operations Research, 54:150–168.

112
Bienstock, D. and Ozbay, D. (2008). Computing robust basestock levels. Discrete Opti-
mization, 5:389–414.

Bookbinder, J. H. and Tan, J.-Y. (1998). Strategies for the probabilistic lot-sizing problem
with service-level constraints. Management Science, 34:1096–1108.

Boute, R. N., Disney, S. M., Lambrecht, M. R., and Van Houdt, B. (2014). Coordinat-
ing lead times and safety stocks under aurocorrelated demand. European Journal of
Production Research, 232:52–63.

Brahimi, N., Absi, N., Dauzère-Pérès, S., and Nordli, A. (2017). Single-item dynamic
lot-sizing problems: An updated survey. European Journal of Operational Research,
263(3):838–863.

Brahimi, N., Dauzere-Peres, S., Najid, N. M., and Nordli, A. (2006). Single item lot sizing
problems. European Journal of Operational Research, 168(1):1–16.

Brandimarte, P. (2006). Multi-item capacitated lot-sizing with demand uncertainty. In-

ternational Journal of Production Research, 44(15):2997–3022.

Brockwell, P. J. and Davis, R. A. (2016). Introduction to Time Series and Forecasting.

Springer Texts in Statistics. Springer International Publishing, Cham.

Campbell, T. and How, J. P. (2015). Bayesian nonparametric set construction for robust
optimization. In 2015 American Control Conference (ACC), pages 4216–4221.

Charnes, J. M., Marmorstein, H., and Zinn, W. (1995). Safety stock determination with
serially correlated demand in a periodic-review inventory system. The Journal of the
Operational Research Society, 46(8):1006–1013.

Chen, X., Sim, M., and Sun, P. (2007). A robust optimization perspective on stochastic
programming. Operation Research, 55:1058–1071.

Christopher, M. (2005). Logistics and Supply Chain Management: Creating Value-Adding

Networks (3rd ed.). FT Prentice Hall.

Coniglio, S., Koster, A. M. C. A., and Spiekermann, N. (2018). Lot sizing with storage
losses under demand uncertainty. Journal of Combinatorial Optimization, 36:763–788.

113
Daneshvari, H. and Shafaei, R. (2021). A new correlated polyhedral uncertainty set for
robust optimization. Computers & Industrial Engineering, 157:107346.

Daryalal, M., Arslan, A. N., and Bodur, M. (2023). Two-stage and lagrangian dual
decision rules for multistage adaptive robust optimization. hal-04090602v1. 09/19/2023.

Disney, S., Farasyn, I., Lambrecht, M., Towill, D., and de Velde, W. V. (2006). Taming
the bullwhip effect whilst watching customer service in a single supply chain echelon.
European Journal of Operational Research, 173(1):151–172.

Dolgui, A. and Proth, J.-M. (2010). Supply Chain Engineering: Useful Methods and
Techniques. Springer London.

Duc, T. T. H., Luong, H. T., and Kim, Y.-D. (2008). A measure of bullwhip effect in
supply chains with a mixed autoregressive-moving average demand process. European
Journal of Operational Research, 187(1):243–256.

El Ghaoui, L., Oustry, F., and Lebret, H. (1998). Robust solutions to uncertain semidef-
inite programs. SIAM Journal on Optimization, 9(1):33–52.

Eppen, G. D. and Martin, R. K. (1988). Determining safety stock in the presence of

stochastic lead time and demand. Management Science, 34(11):1380–1390.

Erkip, N., Hausman, W. H., and Nahmias, S. (1990). Optimal centralized ordering poli-
cies in multi-echelon inventory systems with correlated demand. Management Science,
36(3):381–392.

Fischetti, M. and Monaci, M. (2012). Cutting plane versus compact formulations for un-
certain (integer) linear programs. Mathematical Programming Computation, 4:239–273.

Gabrel, V., Murat, C., and Thiele, A. (2014). Recent advances in robust optimization:
An overview. European Journal of Operational Research, 235(3):471–483.

Gallego, G. (1992). A minmax distribution free procedure for the (q, r) inventory model.
Operations Research Letters, 11(1):55–60.

Gicquel, C., Minoux, M., and Dallery, Y. (2008). Capacitated lot sizing models: a litera-
ture review. hal-00255830.

114
Goerigk, M., Deghdak, K., and T’Kindt, V. (2015). A two-stage robustness approach
to evacuation planning with buses. Transportation Research Part B: Methodological,
78:66–82.

Goerigk, M. and Kurtz, J. (2023). Data-driven robust optimization using deep neural
networks. Computers & Operations Research, 151:106087.

Gumte, K. M., Devi Pantula, P., Miriyala, S. S., and Mitra, K. (2021a). Achieving wealth
from bio-waste in a nationwide supply chain setup under uncertain environment through
data driven robust optimization approach. Journal of Cleaner Production, 291:125702.

Gumte, K. M., Devi Pantula, P., Miriyala, S. S., and Mitra, K. (2021b). Data driven
robust optimization for handling uncertainty in supply chain planning models. Chemical
Engineering Science, 246:116889.

Hamed Mamani, Shima Nassiri, M. R. W. (2016). Closed-form solutions for robust in-
ventory management. Management Science, 63(5):1625–1643.

Han, B., Shang, C., and Huang, D. (2021). Multiple kernel learning-aided robust op-
timization: Learning algorithm, computational tractability, and usage in multi-stage
decision-making. European Journal of Operational Research, 292:1004–1018.

Harris, F. W. (1913). How many parts to make at once. Factory: The Magazine of
Manage- ment, 10(2):135–136.

Huang, K. and Küçüyavuz, S. (2008). On stochastic lot-sizing problems with random lead
times. Operations Research Letters, 36(3):303–308.

Jalilvand-Nejad, A., Shafaei, R., and Shahriari, H. (2016). Robust optimization under
correlated polyhedral uncertainty set. Computers & Industrial Engineering, 92:82–94.

José Alem, D. and Morabito, R. (2012). Production planning in furniture settings via
robust optimization. Computers & Operations Research, 39(2):139–150.

Karimi, B. and Fatemi, G. (2003). The capacitated lot sizing problem: a review of models
and algorithms. Omega, 31:365–378.

Kaufman, L. and Rousseeuw, P. J. (1990). Partitioning Around Medoids (Program PAM).

In Finding Groups in Data, pages 68–125. John Wiley & Sons, Ltd.

115
Keogh, E. and Mueen, A. (2017). Curse of Dimensionality. In Sammut, C. and Webb,
G. I., editors, Encyclopedia of Machine Learning and Data Mining, pages 314–315.
Springer US, Boston, MA.

Lin, X., Zhao, L., Shang, C., He, W., Du, W., and Qian, F. (2022). Data-driven robust
optimization for cyclic scheduling of ethylene cracking furnace system under uncertainty
based on kernel learning. Chemical Engineering Science, 260.

Loger, B., Dolgui, A., Lehuédé, F., and Massonnet, G. (2022). Improving the tractability
of svc-based robust optimization. IFAC-PapersOnLine, 55(10):719–724.

Luong, H. T. (2007). Measure of bullwhip effect in supply chains with autoregressive

demand process. European Journal of Operational Research, 180(3):1086–1097.

Melamed, M., Ben-Tal, A., and Golany, B. (2016). On the average performance of the
adjustable ro and its use as an offline tool for multi-period production planning under
uncertainty. Computational Management Science, 13:293–315.

Metzker Soares, P. (2022). L’optimisation robuste pour des problèmes de lotissement dans
un contexte de rendement incertain. PhD thesis, IMT Atlantique.

Metzker Soares, P., Thevenin, S., Adulyasak, Y., and Dolgui, A. (2022). Robust opti-
mization for lot sizing under yield uncertainty. Computers and Operations Research,
149:106025.

Metzker Soares, P., Thevenin, S., Adulyasak, Y., and Dolgui, A. (2023). Adaptive robust
optimization for lot-sizing under yield uncertainty. European Journal of Operational
Research.

Mohseni, S. and Pishvaee, M. S. (2020). Data-driven robust optimization for wastew-

ater sludge-to-biodiesel supply chain design. Computers & Industrial Engineering,
139:105944.

Montemanni, R. (2006). A benders decomposition approach for the robust spanning tree
problem with interval data. European Journal of Operational Research, 174(3):1479–
1490.

Moon, I. and Gallego, G. (1994). Distribution free procedures for some inventory models.
Journal of the Operational Research Society, 45:651–658.

116
Mutapcic, A. and Boyd, S. P. (2009). Cutting-set methods for robust convex optimization
with pessimizing oracles. Optimization Methods and Software, 24:381 – 406.

Mínguez, R. and Casero-Alonso, V. (2019). On the convergence of cutting-plane methods

for robust optimization with ellipsoidal uncertainty sets.

Nabil Absi, S. K.-S. and Dauzère-Pérès, S. (2011). Uncapacitated lot-sizing problem with
production time windows, early productions, backlogs and lost sales. International
Journal of Production Research, 49(9):2551–2566.

Ning, C. and You, F. (2017). Data-driven adaptive nested robust optimization: General
modeling framework and efficient computational algorithm for decision making under
uncertainty. AIChE Journal, 63(9):3790–3817.

Ning, C. and You, F. (2018a). Data-driven adaptive robust optimization framework based
on principal component analysis. In 2018 Annual American Control Conference (ACC),
pages 3020–3025.

Ning, C. and You, F. (2018b). Data-driven decision making under uncertainty integrating
robust optimization with principal component analysis and kernel smoothing methods.
Computers & Chemical Engineering, 112:190–210.

Pachamanova, D. A. (2002). A robust optimization approach to finance. PhD thesis,

Massachusetts Institute of Technology.

Pahl, J., Voss, S., and Woodruff, D. (2007). Production planning with load dependent
lead times: An update of research. Annals of Operations Research, 153:297–345.

Pahl, J., Voß, S., and Woodruff, D. L. (2005). Load Dependent Lead Times — From
Empirical Evidence to Mathematical Modeling, pages 539–554. Physica-Verlag HD.

Pochet, Y. and Wolsey, L. A. (2006). Production Planning by Mixed Integer Programming.

Springer New York.

Poss, M. (2013). Robust combinatorial optimization with variable budgeted uncertainty.

4OR, 11:75–92.

Pätzold, J. and Schöbel, A. (2020). Approximate cutting plane approaches for exact
solutions to robust optimization problems. European Journal of Operational Research,
284:20–30.

117
Quezada, F., Gicquel, C., and Kedad-Sidhoum, S. (2019). Stochastic dual dynamic integer
programming for a multi-echelon lot-sizing problem with remanufacturing and lost sales.
In IEEE International Conference on Control, Decision and Information Technologies
CODIT 2019, pages 1254–1259. IEEE.

Santos, M. C., Poss, M., and Nace, D. (2018). A perfect information lower bound for
robust lot-sizing problems. Annals of Operations Research, 271:887–913.

Scarf, H. (1958). A min-max solution of an inventory problem. Studies in the mathematical

theory of inventory and production, Standford University Press, pages 201–209.

Schölkopf, B., Platt, J. C., Shawe-Taylor, J., Smola, A. J., and Williamson, R. C.
(2001). Estimating the support of a high-dimensional distribution. Neural computation,
13(7):1443–1471.

Schubert, E. and Rousseeuw, P. J. (2019). Faster k-Medoids Clustering: Improving the

PAM, CLARA, and CLARANS Algorithms. In Amato, G., Gennaro, C., Oria, V., and
Radovanović, M., editors, Similarity Search and Applications, pages 171–187, Cham.
Springer International Publishing.

See, C.-T. and Sim, M. (2010). Robust approximation to multiperiod inventory manage-
ment. Operations Research, 58:583–594.

Shang, C., Chen, W.-H., Stroock, A. D., and You, F. (2020). Robust model predictive
control of irrigation systems with active uncertainty learning and data analytics. IEEE
Transactions on Control Systems Technology, 28(4):1493–1504.

Shang, C., Huang, X., and You, F. (2017). Data-driven robust optimization based on
kernel learning. Computers & Chemical Engineering, 106:464–479.

Shang, C. and You, F. (2018). Robust optimization in high-dimensional data space with
support vector clustering. IFAC-PapersOnLine, 51:19–24.

Shang, C. and You, F. (2019). A data-driven robust optimization approach to scenario-

based stochastic model predictive control. Journal of Process Control, 75:24–39.

Shen, F., Zhao, L., Du, W., Zhong, W., and Qian, F. (2020). Large-scale industrial energy
systems optimization under uncertainty: A data-driven robust optimization approach.
Applied Energy, 259:114199.

118
Simchi-Levi, D. (2008). Designing and Managing the Supply Chain: Concepts, Strategies,
and Case Studies. McGraw-Hill/Irwin.

Solyalı, O., Cordeau, J.-F., and Laporte, G. (2016). The impact of modeling on robust
inventory management under demand uncertainty. Management Science, 62(4):1188–
1201.

Song, J.-S. and Zipkin, P. (2003). Supply chain operations: Assemble-to-order systems.
Handbooks in Operations Research and Management Science, 11:561–596.

Soyster, A. L. (1973). Technical note — convex programming with set-inclusive constraints

and applications to inexact linear programming. Operations Research, 21(5):1154–1157.

Stadtler, H., Kilger, C., and Meyr, H. (2015). Supply Chain Management and Advanced
Planning. Springer Berlin, Heidelberg.

Syntetos, A., Boylan, J., and Disney, S. (2008). Journal of operations research society.
50th Annual Conference of the Operational Research Society 2008, OR50, 60:149–160.

Sözüer, S. and Thiele, A. C. (2016). The State of Robust Optimization. In Robust-

ness Analysis in Decision Aiding, Optimization, and Analytics, International Series
in Operations Research & Management Science, pages 89–112. Springer International
Publishing.

Taş, D., Gendreau, M., Jabali, O., and Jans, R. (2019). A capacitated lot sizing problem
with stochastic setup times and overtime. European Journal of Operational Research,
273(1):146–159.

Thevenin, S., Ben-Ammar, O., and Brahimi, N. (2022). Robust optimization approaches
for purchase planning with supplier selection under lead time uncertainty. European
Journal of Operational Research, 303(3):1199–1215.

Thorsen, A. and Yao, T. (2017). Robust inventory control under demand and lead time
uncertainty. Annals of Operations Research, 257:207–236.

Udom, P. and Phumchusri, N. (2014). A comparison study between time series model
and arima model for sales forecasting of distributor in plastic industry. IOSR Journal
of Engineering, 4:32–38.

119
Varas, M., Maturana, S., Pascual, R., Vargas, I., and Vera, J. (2014). Scheduling produc-
tion for a sawmill: A robust optimization approach. International Journal of Production
Economics, 150:37–51.

Wagner, H. M. and Whitin, T. M. (1958). Dynamic version of the economic lot size model.
Management Science, 5(1):89–96.

Wu, W., Liu, R., Yang, Q., and Quek, T. Q. (2022). Learning-based robust resource
allocation for d2d underlaying cellular network. IEEE Transactions on Wireless Com-
munications, 21(8):6731–6745.

Zeng, B. (2011). Solving Two-stage Robust Optimization Problems by A Constraint-and-

Column Generation Method. Optimization Online.

Zhang, C., Wang, Z., and Wang, X. (2022). Machine learning-based data-driven robust
optimization approach under uncertainty. Journal of Process Control, 115:1–11.

Zhang, Y., Feng, Y., and Rong, G. (2018a). Data-driven rolling-horizon robust opti-
mization for petrochemical scheduling using probability density contours. Computers
& Chemical Engineering, 115:342–360.

Zhang, Y., Jin, X., Feng, Y., and Rong, G. (2018b). Data-driven robust optimization
under correlated uncertainty: A case study of production scheduling in ethylene plant.
Computers & Chemical Engineering, 109:48–67.

Zhao, L. and You, F. (2019). A data-driven approach for industrial utility systems opti-
mization under uncertainty. Energy, 182:559–569.

120
 P ERSONAL PUBLICATIONS LIST

Submitted to an international journal

— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet. Approximate Kernel Learning
Uncertainty Set for Robust Combinatorial Optimization. In press in Informs Journal
On Computing.

International conference with proceedings

— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet. A robust data-driven approach
to supply planning. APMS 2021: IFIP International Conference on Advances in
Production Management Systems, IMT Atlantique, Sep 2021, Nantes, France.

— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet. Improving the tractability of SVC-

based Robust Optimization. MIM 2022: 10th IFAC Conference on Manufacturing
Modelling, Management and Control, Jun 2022, Nantes, France.

International conference without proceedings

— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet. Data-driven robust optimization
approach to supply planning using an approximated uncertainty set. IWLS 2022:
International Workshop on Lot-Sizing, Aug 2022, Oslo, Norway.

National conferences
— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet. Approximation d’un ensemble
d’incertitude pour l’optimisation robuste dirigée par les données. ROADEF 2022
: 23 ème congrès annuel de la Société Française de Recherche Opérationnelle et
d’Aide à la Décision. Villeurbanne - Lyon, France.

121
— B. Loger, A. Dolgui, F. Lehuédé, G. Massonnet, S. Minner, B. Sel. Robust inventory
management under joint demand and lead-time uncertainty. ROADEF 2023 : 24
ème congrès annuel de la Société Française de Recherche Opérationnelle et d’Aide
à la Décision. Rennes, France.

122
Titre : Modèles d’optimisation basés sur les données pour la planification des opérations dans
les Supply Chain industrielles

Mot clés : Optimisation dirigées par les données, Gestion de la Supply Chain, Optimisation
Robuste

Résumé : Face à la complexité croissante des fis en développant différentes méthodes d’op-
chaînes logistiques, l’utilisation d’outils d’aide timisation mathématiques reposant sur l’uti-
à la décision automatisés devient nécessaire lisation de données historiques, dans le but
pour appréhender les multiples sources d’in- de proposer des solutions robustes à plu-
certitude susceptibles de les impacter tout en sieurs problèmes d’approvisionnement et de
garantissant un niveau de performance élevé. planification de production. Des expérimenta-
Pour répondre à ces objectif, les managers ont tions numériques sur des applications de na-
de plus en plus recours à des approches ca- tures diverses permettent de comparer ces
pables d’améliorer la résilience des chaînes nouvelles techniques à plusieurs autres ap-
logistiques en proposant des solutions ro- proches classiques de la littérature pour va-
bustes face aux aléas, pour garantir à la fois la lider leur pertinence en pratique. Les résultats
qualité de service et la maîtrise des coûts liés obtenus démontrent l’intérêt de ces contribu-
à la production, au stockage et au transport tions, qui offrent des performances moyennes
de biens. Alors que la collecte et l’analyse des comparables tout en réduisant leur variabilité
données occupe une place croissante dans la en contexte incertain. En particulier, les so-
stratégie des entreprises, la bonne exploita- lutions restent satisfaisantes lorsqu’elles sont
tion de ces informations afin de caractériser confrontées à des scenarios extrêmes, dont
plus précisément ces incertitudes et leur im- la probabilité d’apparition est faible. Enfin, les
pact sur les opérations devient un enjeu ma- temps de résolution des procédures dévelop-
jeur pour optimiser les systèmes de produc- pées restent compétitifs et laissent envisager
tion et de distribution modernes. Cette thèse la possibilité d’une mise en œuvre sur des cas
se positionne au cœur de ces nouveaux dé- d’application à l’échelle industrielle.
Title: Data driven optimisation models for operations planning in industrial Supply chains

Keywords: Data Driven Optimization, Supply Chain Management, Robust Optimization

Abstract: With the increasing complexity new challenges by developing different math-
of supply chains, automated decision-support ematical optimization methods based on his-
tools become necessary in order to appre- torical data, with the aim of proposing ro-
hend the multiple sources of uncertainty that bust solutions to several supply and produc-
may impact them, while maintaining a high tion planning problems. To validate the prac-
level of performance. To meet these objec- tical relevance of these new techniques, nu-
tives, managers rely more and more on ap- merical experiments on various applications
proaches capable of improving the resilience compare them with several other classical ap-
of supply chains by proposing robust solutions proaches from the literature. The results ob-
that remain valid despite uncertainty, to guar- tained demonstrate the value of these contri-
antee both a quality of service and a con- butions, which offer comparable average per-
trol of the costs induced by the production, formance while reducing their variability in an
storage and transportation of goods. As data uncertain context. In particular, the solutions
collection and analysis become central to de- remain satisfactory when confronted with ex-
fine the strategy of companies, a proper us- treme scenarios, whose probability of occur-
age of this information to characterize more rence is low. Finally, the computational time
precisely these uncertainties and their impact of the procedures developed remain competi-
on operations is becoming a major challenge tive, making them suitable for industrial-scale
for optimizing modern production and distri- applications.
bution systems. This thesis addresses these