Network Flow Models for Robust Binary Optimization with Selective Adaptability

1. Introduction

Robust optimization (RO) has become a well-established paradigm for modeling and solving decision-making problems under uncertainty. It has found diverse applications, particularly in problem settings where the distribution of model parameters is difficult to estimate or when it is necessary to consider worst-case outcomes. Early research focused largely on static RO models, where a single robust solution is applied to all possible parameter realizations (Bertsimas et al. 2013). More recently, there has been significant interest in adaptive RO models, where recourse decisions can be made once additional information is revealed (Yanıkoğlu et al. 2019).

Adaptive RO models, particularly those with discrete recourse variables, are relevant in numerous application domains, such as routing (Eufinger et al. 2020), scheduling (Yan and Kung 2018), facility location (Hanasusanto et al. 2015), network design (Álvarez-Miranda et al. 2015), batching (Bayram et al. 2022), and matching problems (McElfresh et al. 2019). Despite their importance, these models can be significantly more challenging to solve relative to their static counterparts. Specifically, although there are standard techniques to reformulate and solve static RO models as tractable mixed-integer linear programming (MILP) problems, these techniques do not extend directly to adaptive RO models, where discrete recourse decisions must be defined for each and every parameter value lying within an uncertainty set. Few exact solution methods exist for solving these models, and most require specialized and carefully tuned algorithms (Zeng and Zhao 2013, Kämmerling and Kurtz 2020, Arslan and Detienne 2022). Instead, significant attention has been placed on approximation techniques, particularly ones that yield MILP formulations that can be solved directly using standard solvers (see, e.g., Hanasusanto et al. 2015).

In this paper, we consider adaptive robust binary optimization (ARBO) problems with objective uncertainty. We focus, without loss of generality, on ARBO problems with selective adaptability, a term we introduce to describe a class of linking constraints in which first-stage variables fix the values of one or several recourse variables but leave the others unrestricted, allowing the remaining recourse variables to adapt to new information. We show that this structure arises naturally in various planning and sequential decision-making tasks and can also be induced by adding auxiliary variables to general ARBO formulations that do not inherently exhibit this property. More importantly, we demonstrate that the structure of these linking constraints enables the development of new reformulation techniques and versatile solution methods for ARBO problems. Below, we provide a concise summary of our main contributions:

1. We introduce ARBO problems with selective adaptability and show how the structure of the linking constraints can be exploited when we have a description of the convex hull of the recourse feasible space (Section 3). We then leverage ideas from the decision diagram community to obtain such a description, demonstrating that ARBO problems can be reformulated as constrained network flow models (Section 4). In these models, recourse decisions are represented by continuous flows on certain links in a large, capacitated network and where capacity constraints are given by the values of first-stage decisions. Because these network flow models are MILPs, they can be solved directly using standard off-the-shelf solvers.

2. We present new techniques to generate more compact and tractable network flow models to approximate large ARBO problems (Section 5). The first technique utilizes approximate decision diagrams to formulate inner and outer approximations of the recourse feasible space, which enables the generation of both feasible solutions and primal and dual bounds on a given ARBO problem. The second technique, which we term as a generalized multinetwork flow model, pools together a collection of multiple constraints and networks to generate feasible solutions and dual bounds. We show the versatility of these approximation methods and outline new approaches for specifying the size and quality of these approximations.

3. We examine the performance of our exact and approximate network flow models in two sets of numerical experiments (Section 6). For smaller ARBO problems, we find that the exact network flow models can be solved efficiently. For larger ARBO problems, we show that approximate network flow models are essential and can reliably produce near-optimal solutions and high-quality dual bounds in solution times that are orders of magnitude faster than exact approaches and benchmark models. Furthermore, we show that the size and solution times of the approximate models can be reduced substantially while sacrificing very little in terms of solution quality, making these methods highly scalable to large ARBO problems.

All proofs can be found in the Electronic Companion. All vectors and matrices are presented in bold letters, whereas sets are denoted using calligraphic letters. Let $\mathcal{S}$ be a feasible set described using linear and integrality constraints. We use $\text{Relax}(\mathcal{S})$ to describe the relaxation of $\mathcal{S}$ obtained by removing integrality constraints and $\text{Conv}(\mathcal{S})$ to denote the convex hull of $\mathcal{S}$; by definition, $\text{Conv}(\mathcal{S})\subseteq \text{Relax}(\mathcal{S})$. The set of extreme points of any polyhedral set $\mathcal{S}$ is denoted by $\text{Ext}(\mathcal{S})$.

2. Literature Review

We first review the relevant literature on robust binary optimization with objective uncertainty, followed by the literature on decision diagrams (DDs) and network flow models. We refer to Bertsimas et al. (2013) and Yanıkoğlu et al. (2019) for a broader overview of robust optimization.

2.1. Adaptive Robust Binary Optimization

Static robust binary optimization with objective uncertainty is well studied, where robust formulations of polynomially solvable problems can be NP-hard (see, e.g., Buchheim and Kurtz 2018). Nonetheless, these RO problems admit straightforward MILP formulations by reformulating the inner adversarial problem using standard optimality conditions.

In contrast, few algorithms exist for solving ARBO problems with objective uncertainty, where optimal recourse decisions must be defined for any parameter drawn from the uncertainty set. The first and most well-known solution approach is the nested constraint-and-column generation algorithm, which iterates between a master problem and a subproblem that adds new variables and constraints into the master problem (Zeng and Zhao 2013). The second approach is a branch-and-price algorithm where, within each node of a branch-and-bound tree, the algorithm solves a series of iterative pricing problems to add columns that correspond to feasible recourse decisions under a fixed first-stage solution (Arslan and Detienne 2022). Finally, Kämmerling and Kurtz (2020) defined a specialized branch-and-bound algorithm that branches over first-stage solutions while using another algorithm to iteratively refine the dual bound in each node of the tree. The constraint-and-column generation approach is known to scale poorly to larger instances because each iteration of the algorithm adds a copy of all recourse variables and constraints to the master problem (see, e.g., Dumouchelle et al. 2023). On the other hand, the latter two approaches require carefully tuned intermediary steps and may be confined to specific software; for example, many commercial solvers limit the degree of user customization of specific branch-and-bound processes.

The challenges associated with these exact solution methods have motivated considerable interest in approximation techniques (Vayanos et al. 2011, Postek and Hertog 2016). For example, Bertsimas and Georghiou (2015) and Bertsimas and Dunning (2016) considered the use of piecewise constant decision rules to represent the binary recourse variables. Alternatively, Hanasusanto et al. (2015) and Subramanyam et al. (2020) proposed the K-adaptability method, which pre-identifies K recourse solutions to be implemented in the second stage, resulting in an MILP model that has become a commonly used benchmark (see, e.g., Dumouchelle et al. 2023 and Arslan and Detienne 2022).

In this paper, we present a new family of exact and approximate MILP reformulations for ARBO problems with objective uncertainty. We show that these formulations are highly versatile and can produce near-optimal solutions and bounds in significantly less time than existing methods.

2.2. Decision Diagrams and Network Flow Models

Decision diagrams are directed acyclic graphical structures that can be used to model and approximate the feasible space of discrete optimization problems. In recent years, DDs have been used to analyze, solve, and generate bounds and/or cuts for various discrete optimization problems (Bergman et al. 2016). Whereas most of the literature has focused on binary and integer programs, recent work has extended the use of DDs to nonlinear and mixed integer programs (see, e.g., Davarnia and Van Hoeve 2021 and Salemi and Davarnia 2023). With few exceptions, previous work focuses mainly on deterministic problems. Below, we highlight relevant literature on (approximate) decision diagrams and their applications in optimization problems under uncertainty and refer to Castro et al. (2022) for a more comprehensive survey of decision diagrams.

Although DDs can technically be used to model the feasible space of any discrete optimization problem exactly, one limitation is that the size of these diagrams may increase exponentially with problem size. In contrast, approximate DDs can encode inner and outer approximations of the feasible space, where the size of these DDs can be user specified. Many papers examine the tradeoff between the size of the DDs and the quality of the bounds and/or cuts that they generate, especially when compared with standard linear relaxations of various integer programs (Cire and Van Hoeve 2013; Bergman et al. 2014, 2016; Tjandraatmadja and van Hoeve 2019; Serra and Hooker 2020; Davarnia and Van Hoeve 2021). In this context, DDs are typically considered as a way to compete against, or enhance, standard MILP solution methods for deterministic problems and are used in their graphical form with shortest path algorithms. In contrast, we utilize DDs for their network flow representations to generate MILP-based reformulations and approximations of the ARBO problem. Compared with existing solution methods for ARBO problems, these MILP formulations are easier to formulate and solve using standard solvers.

Several works have examined the use of decision diagrams in two-stage stochastic programming with binary recourse decisions (Lozano and Smith 2018, Guo et al. 2021, MacNeil and Bodur 2023). Specifically, a decision diagram is associated with each scenario, making the two-stage problem amenable to standard Benders decomposition techniques that typically apply only to problems with continuous recourse; in the decomposition, shortest path algorithms are applied on the diagrams to generate cuts for the master problem. On the other hand, Serra et al. (2019) considered a monolithic formulation of a two-stage scheduling problem where a subset of constraints are replaced with DD-based network flow constraints. Integrality of the flows across these networks are then enforced using binary variables. Salemi and Davarnia (2023) extended DDs to represent mixed-integer optimization problems, which they used to reformulate the master problem of a Benders decomposition applied to solve the stochastic unit commitment problem. Hooker (2022) proposed the concept of stochastic decision diagrams for stochastic dynamic programming problems, where transition probabilities are associated with arcs of the decision diagram. To our knowledge, Lozano et al. (2022) has been the only application of DDs for robust optimization, where static RO models are reformulated into DD-based constrained shortest path problems for which specialized algorithms can be employed. In contrast to these works, we consider the first use of DD-based network flow models for adaptive robust optimization, which facilitates MILP reformulations and approximations for these problems.

3. Robust Binary Optimization with Selective Adaptability

In this section, we define the ARBO problem of interest, introduce the concept of selective adaptability, and derive several preliminary insights based on the structure of the ARBO problem.

3.1. Problem Definition

Consider an ARBO problem of the form

where $\mathcal{X}\subseteq {\{0,1\}}^m$ defines the feasible set of binary first-stage decisions, $\mathrm{\Xi }≔\{\mathit{\xi }|\mathit{T}\mathit{\xi }\le \mathbf{d}\}$ is a bounded polyhedral uncertainty set, and $\mathcal{Y}\cap \mathcal{S}(\mathbf{x})\subseteq {\{0,1\}}^n$ defines the feasible set of binary recourse decisions where the set $\mathcal{Y}≔\{\mathbf{y}\in {\{0,1\}}^n|\mathit{G}\mathbf{y}\ge \mathbf{h}\}$ captures the constraints that do not depend on $\mathbf{x}$, whereas $\mathcal{S}(\mathbf{x})\subseteq {\mathbb{R}}^n$ models the linking constraints. In this paper, we focus on problems where $\mathcal{S}(\mathbf{x})$ satisfies a condition that we term selective adaptability.

Definition 1

(Selective Adaptability). A linking constraint between first-stage and second-stage decision variables is defined to be selectively adaptive if it can be expressed as $y_i\le x_j$, $y_i=x_j$, or $y_i\ge x_j$ for some $i\in \{1,\dots ,n\}$ and $j\in \{1,\dots ,m\}$. The ARBO problem (1) is said to have selective adaptability if all constraints defining $\mathcal{S}(\mathbf{x})$ are selectively adaptive, that is,

$$\mathcal{S}(\mathbf{x})=\left\{\mathbf{y}\in {\mathbb{R}}^n\hspace{0.17em}|\begin{array}{ll}{y}_i\le x_j & \forall (i,j)\in {\mathcal{U}}_1\\ y_i=x_j& \forall (i,j)\in {\mathcal{U}}_2\\ y_i\ge x_j& \forall (i,j)\in {\mathcal{U}}_3\end{array}\right\},$$(2)

where ${\mathcal{U}}_1,\hspace{0.17em}{\mathcal{U}}_2,\hspace{0.17em}{\mathcal{U}}_3$ are disjoint sets, and ${\mathcal{U}}_1\cup {\mathcal{U}}_2\cup {\mathcal{U}}_3\subseteq \{1,\dots ,n\}\times \{1,\dots ,m\}$.

An ARBO problem with selective adaptability implies that the value of each binary recourse variable $y_i\in \{0,1\}$ is either (i) fixed by the values of first-stage decisions $\mathbf{x}\in \mathcal{X}$ (because of constraints corresponding to ${\mathcal{U}}_2$ and those associated with ${\mathcal{U}}_1$ and ${\mathcal{U}}_3$ if $x_j$ is 0 or 1, respectively) or (ii) not restricted by the values of first-stage decisions $\mathbf{x}\in \mathcal{X}$ (because the constraints associated with $y_i$ become redundant when ${\mathcal{U}}_2=\varnothing$ and when $x_j=1$ for all $(i,j)\in {\mathcal{U}}_1$ and $x_j=0$ for all $(i,j)\in {\mathcal{U}}_3$). Practically, this implies that individual recourse variables not impacted by the first-stage decisions can be used to “adapt” to new information.

3.2. Models with Selective Adaptability

We give a few examples of ARBO problems with selective adaptability and then show that our examination of this problem structure comes without loss of generality. We begin with two examples of sequential decision-making problems that naturally exhibit this linking constraint structure:

• In a facility location problem, a planner can make an assignment (or shipping) decision between a potential location j and an individual customer l only if a facility (e.g., hub, warehouse, factory) has been placed at location j (Kämmerling and Kurtz 2020). Assuming that customer demand is given by $\mathit{\xi }$, the sequential nature of this decision-making process can be represented by linking constraints of the form $y_{lj}\le x_j$, where $y_{lj}$ is a particular location-customer assignment decision and $x_j$ determines whether a facility is placed in location j.

• In a project investment problem, an investor with a limited budget aims to maximize investment returns across two decision stages by investing in different projects where early-stage investments come with risk but higher return rates (Arslan and Detienne 2022). For each project $i\in \{1,\dots ,n\}$, the linking constraint $y_i\ge x_i$ captures early-stage commitment decisions. We revisit this problem in our numerical experiments in Section 6.

The concept of selective adaptability is also relevant when relaxing static robust optimization problems. Specifically, note that a static RO problem can be written as

where $\mathcal{X}=\mathcal{Y}$ defines the set of feasible solutions and $\mathbf{y}=\mathbf{x}$ removes the ability to adapt to any new information. In this context, replacing $\mathbf{y}=\mathbf{x}$ with $\mathbf{y}\le \mathbf{x}$, and relaxing $\mathcal{X}$ so that $\mathcal{X}\supset \mathcal{Y}$, provides a natural framework for introducing some flexibility into static models. When $\mathbf{x}$ and $\mathbf{y}$ are both binary vectors, this adaptive model can be viewed as one that preselects or narrows down the set of decisions to be considered before more information becomes available. In operational contexts, this can be useful for planning, communication, or training purposes. We revisit this in Section 6.

3.3. Preliminary Model Insights

One possible reformulation of Problem (1) is as an infinite-dimensional MILP model, that is,

where v represents the worst-case recourse objective value and where a set of variables (${\mathbf{y}}^{\xi }$) and constraints must be defined for every $\mathit{\xi }\in \mathrm{\Xi }$. Zeng and Zhao (2013) proposed a nested constraint-and-column generation algorithm to iteratively refine an approximation of (3), where the set $\mathrm{\Xi }$ is replaced with a finite set $\widehat{\mathrm{\Xi }}\subset \mathrm{\Xi }$ that is iteratively enlarged through a subproblem. However, given that (3) is an MILP, the approach can quickly become intractable as the size of $\widehat{\mathrm{\Xi }}\subset \mathrm{\Xi }$ is increased.

Arslan and Detienne (2022) observed that Problem (1) permits an alternative reformulation as

This observation can be made by noting that the inner minimization of Problem (1) is equivalent to optimizing over $\text{Conv}(\mathcal{Y}\cap \mathcal{S}(\mathbf{x}))$. Because this set is convex, we can apply the minimax theorem (Neumann 1928) to swap the order of the inner max and min operators, which leads to model (4).

We now show that there exists an even simpler reformulation when $\mathcal{S}(\mathbf{x})$ satisfies selective adaptability. This reformulation is shown in the next proposition, which relies on the following lemma.

Lemma 1.

If Problem (1) has selective adaptability, then

$$\text{Conv}(\mathcal{Y}\cap \mathcal{S}(\mathbf{x}))=\text{Conv}(\mathcal{Y})\cap \mathcal{S}(\mathbf{x}),\hspace{0.17em}\forall \mathbf{x}\in {\{0,1\}}^m.$$

Proposition 1.

If Problem (1) has selective adaptability, then model (3) can be reformulated as shown, where $\mathbf{y}$ becomes a vector of continuous variables:

$$\underset{\mathbf{x},\mathbf{y},\mathit{\lambda }}{\text{min}}\hspace{0.17em}\{{\mathbf{c}}^{\mathrm{\top }}\mathbf{x}+{\mathbf{d}}^{\mathrm{\top }}\mathit{\lambda }|{\mathit{T}}^{\mathrm{\top }}\mathit{\lambda }=\mathbf{y},\hspace{0.17em}\mathbf{y}\in \text{Conv}\left(\mathcal{Y}\right),\hspace{0.17em}\mathbf{y}\in \mathcal{S}\left(\mathbf{x}\right),\hspace{0.17em}\mathbf{x}\in \mathcal{X},\hspace{0.17em}\mathit{\lambda }\ge \mathbf{0}\}.$$(5)

The second and third constraints come directly from applying Lemma 1 to model (4), whereas the term $\max \{{\mathit{\xi }}^{\mathrm{\top }}\mathbf{y}|\mathit{\xi }\in \mathrm{\Xi }\}$ in model (4) is reformulated using strong duality. Namely, recall that because $\mathrm{\Xi }≔\{\mathit{\xi }|\mathit{T}\mathit{\xi }\le \mathbf{d}\}$, the dual of this maximization term can be rewritten as $\min \{{\mathbf{d}}^{\mathrm{\top }}\mathit{\lambda }|{\mathit{T}}^{\mathrm{\top }}\mathit{\lambda }=\mathbf{y},\mathit{\lambda }\ge \mathbf{0}$}. Thus, the new variables $\mathit{\lambda }$ in model (5) correspond to the dual variables of this reformulation.

Proposition 1 is important for two reasons and motivates the rest of this paper. First, model (5) is equivalent to a static RO problem, that is, where both $\mathbf{x}\in {\{0,1\}}^m$ and $\mathbf{y}\in {\mathbb{R}}^n$ represent “first-stage” decision variables. Practically, this implies that the original ARBO problem can be solved as an MILP if we have a polyhedral representation of $\text{Conv}(\mathcal{Y})$. Section 4 focuses on representing $\text{Conv}(\mathcal{Y})$ in an extended network-based space. Second, Proposition 1 provides intuition for developing approximation methods, which are the focus of Section 5. In particular, replacing $\text{Conv}(\mathcal{Y})$ with an inner (or outer) approximation of it will generate a valid primal (or dual) bound on the optimal value of (5). For example, replacing $\text{Conv}(\mathcal{Y})$ with the continuous relaxation of $\mathcal{Y}$, that is,

generates a lower bound on the optimal value of model (5) because $\text{Conv}(\mathcal{Y})\subseteq \text{Relax}(\mathcal{Y})$.

For the rest of this paper, we assume the ARBO problem has selective adaptability. This assumption is made for ease of exposition and comes without loss of generality. In particular, any ARBO problem without selective adaptability can be reformulated into one that has this property by adding auxiliary variables into the recourse problem. Specifically, let $\widehat{\mathbf{F}}\mathbf{x}+\widehat{\mathbf{G}}\mathbf{y}\le \widehat{\mathbf{h}}$ be a set of linking constraints that are not in the form of (2), and let $\widehat{\mathcal{I}}\subseteq \{1,\dots ,m\}$ be the indices of the $\mathbf{x}$ variables that appear in these constraints. Now, let ${\mathit{y}}^{\text{aux}}\in {\{0,1\}}^{|\widehat{\mathcal{I}}|}$ be a new set of auxiliary recourse variables. For each $j\in \widehat{\mathcal{I}}$, we can replace $x_j$ in $\widehat{\mathbf{F}}\mathbf{x}+\widehat{\mathbf{G}}\mathbf{y}\le \widehat{\mathbf{h}}$ with the corresponding $y_j^{\text{aux}}$ variable, making these constraints a part of $\mathcal{Y}$, and redefine the linking constraints as $\mathcal{S}(\mathbf{x})=\{y_j^{\text{aux}}=x_j,\hspace{0.17em}\forall j\in \widehat{\mathcal{I}}\}$. This technique is inspired by Arslan and Detienne (2022) and implies that all results derived in our paper are applicable to any ARBO problem with or without selective adaptability.

We conclude this section with two additional remarks on the generalizability of upcoming results.

Remark 1.

The reformulations and approximations that are proposed in this paper apply to problems where first-stage decisions are mixed-integer, that is, $\mathcal{X}\subseteq {\{0,1\}}^{m_1}\times {\mathbb{Z}}^{m_2}\times {\mathbb{R}}^{m_3}$, as long as only the binary variables defining $\mathcal{X}$ appear in the linking constraints $\mathcal{S}(\mathbf{x})$.

Remark 2.

After Section 4, one may wonder whether we can use the network flow formulations to directly represent $\text{Conv}(\mathcal{Y}\cap \mathcal{S}(\mathbf{x}))$ in constraint (4b) by building a decision diagram representation of the set $(\mathbf{x},\mathbf{y})\in \mathcal{Y}\cap \mathcal{S}(\mathbf{x})$ rather than artificially inducing selective adaptability using new auxiliary variables. We note that these two ideas are in fact equivalent. A decision diagram for the decision region $(\mathbf{x},\mathbf{y})\in \mathcal{Y}\cap \mathcal{S}(\mathbf{x})$ must by definition have m more layers than one for $\mathbf{y}\in \mathcal{Y}$ (see Section 4). This would be exactly the same as the decision diagram for the set $({\mathbf{y}}^{\text{aux}},\mathbf{y})\in \mathcal{Y}\cap \mathcal{S}({\mathbf{y}}^{\text{aux}})$, where ${\mathbf{y}}^{\text{aux}}\in {\{0,1\}}^m$ are auxiliary variables introduced to reformulate a general ARBO problem into one with selective adaptability (where $\mathbf{x}={\mathbf{y}}^{\text{aux}}$ become the linking constraints).

4. Constrained Network Flow Reformulations

In this section, we present an MILP formulation for model (4) when $\text{Relax}(\mathcal{Y})\ne \text{Conv}(\mathcal{Y})$. Specifically, in Section 4.1, we introduce DDs and the network flow formulation used to obtain a polyhedral description of $\text{Conv}(\mathcal{Y})$. Then, in Section 4.2, we integrate the formulation with first-stage decisions.

4.1. Reformulating the Recourse Feasible Space

We first describe a general procedure for obtaining a decision diagram encoding of the feasible solutions in $\mathcal{Y}$ and then represent this diagram using a network flow formulation.

4.1.1. Decision Diagram Formulation.

A binary decision diagram (BDD) is a graphical structure that can be used to encode the feasible space of a binary optimization problem (Bergman et al. 2016). Specifically, a BDD is a directed acyclic graph $\mathcal{D}=(\mathcal{V},\mathcal{A})$ with nodes $\mathcal{V}$ and arcs $\mathcal{A}$. The nodes are partitioned into $n+1$ nonempty layers $\mathcal{V}=({\mathcal{V}}_1,\dots ,{\mathcal{V}}_{n+1})$, whereas the directed arcs connect nodes in adjacent layers from ${\mathcal{V}}_i$ to ${\mathcal{V}}_{i+1}$ for $i\in \{1,\dots ,n\}$. The sets ${\mathcal{V}}_1=\{\mathrm{r}\}$ and ${\mathcal{V}}_{n+1}=\{\mathrm{t}\}$ are each composed of a single node, defined as the root node and terminal node, respectively. Each arc in the network has a label of either zero or one, and ${\mathcal{A}}_i^0$ and ${\mathcal{A}}_i^1$ define the subset of zero and one arcs leaving nodes in ${\mathcal{V}}_i$, respectively. A decision diagram $\mathcal{D}$ is a valid representation of a feasible space $\mathcal{Y}$ if each path from root node r to terminal node t in $\mathcal{D}$ can be mapped to a solution $\mathbf{y}\in \mathcal{Y}$ and vice versa. This mapping of paths to solution is defined by the zero-one label on each arc in the path. Specifically, for any arc j in the path, if $j\in {\mathcal{A}}_i^0$, then $y_i=0$, or if $j\in {\mathcal{A}}_i^1$, then $y_i=1$. For example, a decision diagram of $\mathcal{Y}≔\{\mathbf{y}\in {\{0,1\}}^5|y_1+y_2+2y_3+2y_4+3y_5\le 4\}$ is given in Figure 1. Finally, we note that the definition of BDDs in the literature generally includes arc weights that correspond to the value of objective coefficients of an optimization problem. Because these coefficients are not deterministic in our problem, we instead consider only “unweighted” decision diagrams.

DDs are closely related to the state-transition graph in the dynamic programming literature (Hooker 2013). In particular, nodes and arcs in the DD can be mapped to “states” and feasible “actions” of a recursive formulation where decisions are made sequentially. Many binary optimization problem structures admit simple recursive formulations that can be used to obtain decision diagrams (Bergman et al. 2016, 2022). In general, a recursive formulation of a deterministic binary optimization problem can be written as a Bellman equation of the form

where $\mathit{S}$ denotes the state of the system, ${\mathcal{Q}}_i(\mathit{S})\subseteq \{0,1\}$ denotes the set of feasible actions at stage i for state $\mathit{S}$, $T_i(\mathit{S},y_i)$ defines the new state after taking action $y_i$, and the pair $V_i(\mathit{S})$ and $f_i(\mathit{S},y_i)$ are used to capture the long-term values and immediate rewards of taking specific actions, respectively. However, because we need to generate only unweighted DDs, we need only the state-transition graph corresponding to the states $\mathit{S}$ and feasible actions $Q_i(\mathit{S})$. To define the diagram, we first modify $T_n(\mathit{S},y_i)$ to be $T_n(\mathit{S},y_i)=\{\mathrm{t}\}$ and then map each state to a node and each feasible action $y_i\in Q_i(\mathit{S})$ to a unique arc with label $y_i$ that links $\mathit{S}$ and $T_i(\mathit{S},y_i)$. This creates a diagram with a single root node r and terminal node t, where each path from r to t corresponds to a sequence of actions in the recursive formulation.

Example 1.

Consider a feasible set $\mathcal{Y}=\{\mathbf{y}\in {\{0,1\}}^n|{\sum }_{i=1}^n{g}_i{y}_i\le h\}$ that is defined by a knapsack constraint. The recursive formulation of $\mathcal{Y}$ is defined by feasible actions ${\mathcal{Q}}_i(\mathit{S})=\{y_i\in \{0,1\}|\mathit{S}+g_i{y}_i\le h\}$ and state-transition function $T_i(\mathit{S},y_i)=\mathit{S}+g_i{y}_i$ with an initial state ${\mathit{S}}_1=0$. Each stage i in the recursive formulation thus corresponds to a layer i in the diagram. Each state in layer i represents the total amount of capacity used by the selection of items among $1,\dots ,i$, and feasible actions correspond to whether an additional item $i+1$ can be placed in the knapsack given the current state. To obtain a single terminal node in the decision diagram, we merge all nodes in layer $n+1$, which is equivalent to replacing ${\mathrm{T}}_n(\mathit{S},y_n)=\mathit{S}+g_n{y}_n$ with ${\mathcal{T}}_n(\mathit{S},y_n)=h$ in the recursive formulation. Note that this change in the state-transition matrix at the final stage n does not change the set of feasible actions in the recursive formulation.

Finally, we note that DDs can typically be reduced to encode the same set of feasible solutions with far fewer nodes and arcs. Figure 1 is a reduced DD for the knapsack set defined earlier. To obtain reduced DDs, we can use a simple bottom-up merging technique (Bryant 1992, Bergman et al. 2016). Starting with nodes in layer n, we can merge any nodes that have the same set of outgoing arc types and destinations, repeating this procedure from layers $n-1$ to 1.

4.1.2. Recourse Network Flow Formulation.

Let $\mathcal{D}$ denote a DD representation of $\mathcal{Y}$. Following the notation of Castro et al. (2022), we use $\text{NF}(\mathcal{D})$ to denote the network flow model of $\mathcal{D}$, which relates arc flows $\mathbf{z}$ to values of $\mathbf{y}\in \mathcal{Y}$ in the recourse problem, that is,

where $\mathit{A}\in {\{-1,0,1\}}^{|\mathcal{V}|\times |\mathcal{A}|}$ denotes the node-arc incidence matrix and $\mathbf{b}$ a vector of zeros with the exception of $b_1=-1$ and $b_{|\mathcal{V}|}=1$. The constraint set $\mathit{A}\mathbf{z}=\mathbf{b}$ defines standard flow conservation constraints. The second set of constraints link the value of each variable $y_i$ to the sum of total flow over the one-arcs ${\mathcal{A}}_i^1$ in layer $i\in \{1,\dots ,n\}$.

4.2. A Complete Network Flow Reformulation

A key property of $\text{NF}(\mathcal{D})$ is that its projection onto the variables $\mathbf{y}$, denoted by ${\text{Proj}}_{\mathbf{y}}(\text{NF}(\mathcal{D}))$, is equal to the convex hull of $\mathcal{Y}$ (Tjandraatmadja and van Hoeve 2019, Davarnia and Van Hoeve 2021). We can use this property to derive an MILP formulation of Problem (1), as shown next.

Proposition 2.

Problem (1) can be reformulated into the constrained network flow problem

$$\underset{\mathbf{x},\mathbf{y},\mathbf{z},\mathit{\lambda }}{\text{min}}\hspace{0.17em}\{{\mathbf{c}}^{\mathrm{\top }}\mathbf{x}+{\mathbf{d}}^{\mathrm{\top }}\mathit{\lambda }|{\mathit{T}}^{\mathrm{\top }}\mathit{\lambda }=\mathbf{y},\hspace{0.17em}(\mathbf{y},\mathbf{z})\in \text{NF}\left(\mathcal{D}\right),\hspace{0.17em}\mathbf{y}\in \mathcal{S}\left(\mathbf{x}\right),\hspace{0.17em}\mathbf{x}\in \mathcal{X},\hspace{0.17em}\mathit{\lambda }\ge \mathbf{0}\}.$$(8)

Model (8) is an exact MILP reformulation of any ARBO problem with selective adaptability. It is derived by replacing the constraint $\mathbf{y}\in \text{Conv}(\mathcal{Y})$ in model (5) with the network flow formulation in (7). As we show in Section 6, this model can be tractably solved for problems of smaller sizes.

Example 2.

Consider the adaptive robust knapsack problem

$$\underset{\mathbf{x}\in {\{0,1\}}^5}{\text{min}}\hspace{0.17em}\underset{\mathit{\xi }\in \mathrm{\Xi }}{\text{max}}\hspace{0.17em}\underset{\mathbf{y}}{\text{min}}\hspace{0.17em}\{{\mathbf{c}}^{\mathrm{\top }}\mathbf{x}+{\mathit{\xi }}^{\mathrm{\top }}\mathbf{y}|\mathbf{y}\le \mathbf{x},\hspace{0.17em}\mathbf{y}\in \mathcal{Y}≔\{\mathbf{y}\in {\{0,1\}}^5|y_1+y_2+2y_3+2y_4+3y_5\le 4\left\}\right\}.$$(9)

Recall that the DD for $\mathcal{Y}$ is illustrated in Figure 1. Suppose $\mathrm{\Xi }=\left\{\mathit{\xi }|\mathit{\xi }\ge {\mathit{\xi }}^0,\hspace{0.17em}{\Vert \mathit{\xi }-{\mathit{\xi }}^0\Vert }_1\le \delta \right\}$. Then, problem (9) can be reformulated into model (8), where $\mathit{A}$ is the node-arc incidence matrix of the decision diagram given in Figure 1, and

$$\mathit{T}=\left[\begin{array}{ccccc}-1& 0& 0& 0& 0\\ 0& -1& 0& 0& 0\\ 0& 0& -1& 0& 0\\ 0& 0& 0& -1& 0\\ 0& 0& 0& 0& -1\\ 1& 1& 1& 1& 1 \end{array}\right],\hspace{0.17em}\mathbf{b}=\left[\begin{array}{c}-1\\ 0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right],\hspace{0.17em}\mathbf{d}=\left[\begin{array}{c}-{\xi }_1^0\\ -{\xi }_2^0\\ -{\xi }_3^0\\ -{\xi }_4^0\\ -{\xi }_5^0\\ {({\mathit{\xi }}^0)}^{\mathrm{\top }}\mathbf{1}+\delta \end{array}\right],\hspace{1em}\begin{array}{ll}{y}_1=z_2& \\ y_2=z_4+z_6& \\ y_3=z_8+z_{10}+z_{12}& \\ y_4=z_{14}+z_{16}& \\ y_5=z_{19}.& \end{array}$$

The indices $\{2,4,6,8,10,12,14,16,19\}$ correspond to the solid lines in Figure 1 when labeled from left to right in each layer, starting with the first layer.

5. Network Flow Approximations

In this section, we present methods for constructing approximate decision diagrams, which in turn yield more tractable and compact network flow models. We then demonstrate that these models can generate both solutions and bounds for ARBO problems.

5.1. Approximate Decision Diagrams

A restricted decision diagram of $\mathcal{Y}$ contains paths that map to a strict subset of the feasible solutions in $\mathcal{Y}$. On the other hand, a relaxed decision diagram contains paths that map to a superset of solutions that include all solutions in $\mathcal{Y}$. Restricted and relaxed decision diagrams can be generated by merging states in the recursive formulation of $\mathcal{Y}$ (Castro et al. 2022). Next, we outline two general “top-down” approaches for building restricted and relaxed decision diagrams.

5.1.1. Merging with Width-Based Thresholds.

In Algorithm 1 we present a common approach that merges nodes whenever a threshold on the diagram’s “width” has been exceeded. Specifically, each layer in the diagram is constrained to have a width of at most W; that is, there can be at most W nodes per layer. Considering the recursive formulation, this width constraint is equivalent to having at most W number of different states at stage i. Note that if $W=\infty$ in Algorithm 1, then an exact decision diagram will be generated.

Algorithm 1

(Width-Based Decision Diagram Construction Procedure)

Input: A recursive formulation of $\mathcal{Y}$, a width parameter W

Output: A decision diagram $\mathcal{D}$

• 1: Initialization: let ${\mathit{S}}_1(\mathrm{r})$ denote the initial state of the system (at root node r), ${\mathcal{V}}_1=\{\mathtt{r}\},{\mathcal{V}}_2,\dots ,{\mathcal{V}}_n=\varnothing ,\hspace{0.17em}{\mathcal{V}}_{n+1}=\{\mathtt{t}\},\hspace{0.17em}\mathcal{A}=\varnothing$

• 2: for $i\in \{1,\dots ,n-1\}$, do

• 3:  for $u\in {\mathcal{V}}_i$, do

• 4:   for $y_i\in Q_i(\mathit{S}(u))$, do

• 5:    if $\exists \hspace{0.17em}{u}^{\prime }\in {\mathcal{V}}_{i+1}$ where $\mathit{S}(u^{\prime })=T_i(\mathit{S},y_i)$, then add arc from u to $u^{\prime }$ with label $y_i$

• 6:    else add node $u^{\prime }$ with state $T_i(\mathit{S},y_i)$ to ${\mathcal{V}}_{i+1}$ and add arc from u to $u^{\prime }$ with label $y_i$

• 7:  if $|{\mathcal{V}}_{i+1}|>W$, then

• 8:   ${\widehat{\mathcal{V}}}_{i+1}\leftarrow \mathsf{\text{select}}({\mathcal{V}}_{i+1},|{\mathcal{V}}_{i+1}|-W+1)$

• 9:   $\mathsf{\text{merge}}({\widehat{\mathcal{V}}}_{i+1})$

• 10: for $u\in {\mathcal{V}}_n$, do

• 11:   for $y_i\in Q_i(\mathit{S}(u))$, do add arc from u to t with label $y_i$

• 12: Reduce $\mathcal{D}=(\mathcal{V},\mathcal{A})$

• 13: return $\mathcal{D}$

In step 8, the algorithm requires a rule for selecting $(\left|{\mathcal{V}}_{i+1}\right|-W+1)$ nodes from the set ${\mathcal{V}}_{i+1}$. The simplest and most commonly used procedure is random selection (although designing node selection procedures has become a more active topic of research in recent years; see, e.g., van Hoeve 2022). Once this subset of nodes ${\widehat{\mathcal{V}}}_i$ has been selected, they can be merged to create either restricted or relaxed DDs. To create relaxed DDs, merging operators must assign a new state to each merged node such that all subsequent feasible actions in the original recursive formulation remain feasible under this new state (Hooker 2013). Similarly, a restricted DD can be created by assigning a new state to the merged node such that a subset of feasible actions in the original recursive are retained without introducing infeasible actions. We give an example below.

Example 3.

Consider a feasible space defined by a single knapsack constraint (see Example 1). A relaxed (restricted) decision diagram can be generated by merging the nodes ${\widehat{\mathcal{V}}}_{i+1}$ and assigning the new node a state that is the minimum (maximum) value of the states in the merged set.

We note that the node merge operation can be performed while building a layer, rather than after the layer is completely built, to reduce the memory requirement if desired. Furthermore, we note that we can also generate restricted diagrams simply by discarding nodes rather than merging them in Algorithm 1 because we are effectively removing all feasible actions associated with that state-stage combination in the recursive formulation. Finally, step 12 is to reduce the size of the decision diagram using the bottom-up approach described at the end of Section 4.1.1.

5.1.2. Merging with Distance-Based Thresholds.

We now propose a merging approach based on the similarity of states in each layer; that is, we merge two nodes only when their respective state values are within some “distance” of each another. Algorithm 2 outlines this approach, which revolves around a partitioning of nodes in each layer into subgroups, where nodes in each subgroup are then merged. Any suitable set partitioning method can be used, as long as two conditions are met: (i) within each subgroup, some predefined notion of distance between any pair of states does not exceed a user-specified parameter Q, and (ii) for any two subgroups, there exists a pair of nodes, one node in each subgroup, where the distance exceeds Q. The latter condition ensures the fewest number of partitions given Q. To the best of our knowledge, we are the first to outline and implement this distance-based approach for generating approximate decision diagrams.

Algorithm 2

(Distance-Based Decision Diagram Construction Procedure)

Input: A recursive formulation of $\mathcal{Y}$, a distance function $d(·,·)$, and distance parameter Q

Output: A decision diagram $\mathcal{D}$

• 1: Initialization: let ${\mathit{S}}_1(\mathtt{r})$ denote the initial state of the system (at root node r), ${\mathcal{V}}_1=\{\mathrm{r}\},{\mathcal{V}}_2,\dots ,{\mathcal{V}}_n=\varnothing ,\hspace{0.17em}{\mathcal{V}}_{n+1}=\{\mathtt{t}\},\hspace{0.17em}\mathcal{A}=\varnothing$

• 2: for $i\in \{1,\dots ,n-1\}$ do

• 3:  Steps 3–6 from Algorithm 1

• 4:  ${\widehat{\mathcal{V}}}_{i+1,1},{\widehat{\mathcal{V}}}_{i+1,2},\dots \leftarrow \mathsf{\text{partition}}({\mathcal{V}}_{i+1},Q)$

• 5:  for ${\widehat{\mathcal{V}}}_{i+1}\in \{{\widehat{\mathcal{V}}}_{i+1,1},{\widehat{\mathcal{V}}}_{i+1,2},\dots \}$, do

• 6:   $\mathsf{\text{merge}}({\widehat{\mathcal{V}}}_{i+1})$

• 7: Steps 10–13 from Algorithm 1

Example 4.

We give an example of a partition function (step 4 of Algorithm 2). Consider the feasible space defined in Example 3, for which we can define the following node partitioning procedure: (i) order nodes according to state values, from smallest to largest; (ii) starting with the smallest state, create a new subgroup for the next node if and only if its state is more than Q units larger than the smallest state in the current subgroup. Here, distance between states is simply defined as the absolute difference between state values.

The distance-based merging approach is motivated by our initial attempts at using width-based merging in preliminary experiments. We found that it was difficult to finely tune the size of the diagram and quality of the approximation through the use of the width threshold parameter W. In particular, when the width threshold is exceeded, it is often exceeded by a large margin, which typically results in the merging of hundreds or thousands of nodes into a single node. Furthermore, because the width constraint must be satisfied at each layer, there are essentially no restrictions on which nodes can or cannot be merged. In contrast, the key advantage of the distance-based approach is that both the size and the quality of the approximation is closely tied to the value of Q. In distance-based merging, Q determines precisely which nodes in each layer can and cannot be merged based on the similarity of their states. Furthermore, at each layer there may be many subgroups, but the number of nodes to be merged within each subgroup may be small. In general, we find that this distance-based merging procedure leads to approximations of higher quality and more precision in tuning the size and tractability of the corresponding models.

Finally, we remark that the distance-based approach can be customized to generate even more refined approximations by merging subgroups with some probability $p<1$. This additional feature can generate approximations that lie in between those generated solely by incrementally increasing the value of Q (which could be discrete when all states are integer-valued, as in Example 3).

5.2. Generating First-Stage Solutions and Dual Bounds

In this section, we focus on the use of (relaxed) decision diagrams to formulate MILP models that simultaneously generate first-stage solutions and dual bounds for the ARBO problem. A brief discussion of restricted diagrams can be found in Section EC.2 of the Electronic Companion.

Proposition 3.

Let ${\mathcal{D}}_{\mathtt{\text{outer}}}$ denote a relaxed decision diagram of $\mathcal{Y}$. Then, the model

$$\underset{\mathbf{x},\mathbf{y},\mathbf{z},\mathit{\lambda }}{\text{min}}\hspace{0.17em}\{{\mathbf{c}}^{\mathrm{\top }}\mathbf{x}+{\mathbf{d}}^{\mathrm{\top }}\mathit{\lambda }|{\mathit{T}}^{\mathrm{\top }}\mathit{\lambda }=\mathbf{y},\hspace{0.17em}(\mathbf{y},\mathbf{z})\in \text{NF}\left({\mathcal{D}}_{\mathtt{\text{outer}}}\right),\hspace{0.17em}\mathbf{y}\in \text{Relax}\left(\mathcal{Y}\right),\hspace{0.17em}\mathbf{y}\in \mathcal{S}\left(\mathbf{x}\right),\hspace{0.17em}\mathbf{x}\in \mathcal{X}\}$$(10)

generates a feasible solution $\mathbf{x}\in \mathcal{X}$ and a lower bound on the optimal objective value of model (1).

The proposition follows from the observation that ${\mathcal{D}}_{\mathtt{\text{outer}}}$ represents a superset of $\mathcal{Y}$, which implies that $\text{Conv}(\mathcal{Y})\subseteq {\text{Proj}}_{\mathbf{y}}(\text{NF}({\mathcal{D}}_{\mathtt{\text{outer}}}))$. Furthermore, (10) includes the additional term $\mathbf{y}\in \text{Relax}(\mathcal{Y})$ to tighten the dual bound that is generated because, in general, $\text{Relax}(\mathcal{Y})⊈{\text{Proj}}_{\mathbf{y}}(\text{NF}({\mathcal{D}}_{\mathtt{\text{outer}}}))$ and ${\text{Proj}}_{\mathbf{y}}(\text{NF}({\mathcal{D}}_{\mathtt{\text{outer}}}))⊈\text{Relax}(\mathcal{Y})$ (i.e., $\text{Relax}(\mathcal{Y})$ can have fractional extreme points, whereas ${\text{Proj}}_{\mathbf{y}}(\text{NF}({\mathcal{D}}_{\mathtt{\text{outer}}}))$ is an integral polyhedron but includes solutions $\widehat{\mathbf{y}}\in {\{0,1\}}^n$ that are not in $\mathcal{Y}$). The intersection of $\text{Relax}(\mathcal{Y})$ and ${\text{Proj}}_{\mathbf{y}}(\text{NF}({\mathcal{D}}_{\mathtt{\text{outer}}}))$ can thus yield a better outer approximation of $\text{Conv}(\mathcal{Y})$ compared with either sets alone (see Section EC.3.3 for details).

Whereas (10) uses a single decision diagram to generate an outer approximation of $\text{Conv}(\mathcal{Y})$, we propose a generalization using a collection of diagrams and feasible sets.

Proposition 4.

Suppose that $\mathcal{Y}=\{\mathbf{y}\in {\{0,1\}}^n|{({\mathbf{g}}^1)}^{\mathrm{\top }}\mathbf{y}\ge h_1,\dots ,{({\mathbf{g}}^J)}^{\mathrm{\top }}\mathbf{y}\ge h_J\}$, and let $\mathcal{J}=\{1,\dots ,J\}$ denote the set of constraint indices. Now, suppose we are given a collection of subsets of indices ${\mathcal{J}}_1^1,\dots ,{\mathcal{J}}_k^1$ and ${\mathcal{J}}_1^2,\dots ,{\mathcal{J}}_q^2$, where ${\cup }_{i=1}^k{\mathcal{J}}_i^1\subseteq \mathcal{J}$ and ${\cup }_{i=1}^q{\mathcal{J}}_i^2\subseteq \mathcal{J}$. For an arbitrary set ${\mathcal{J}}_i\subseteq \mathcal{J}$, let ${\mathcal{D}}^{{\mathcal{J}}_i}$ and ${\mathcal{D}}_{\mathtt{\text{outer}}}^{{\mathcal{J}}_i}$ denote an exact and relaxed decision diagram for the feasible set described by the constraints with indices in ${\mathcal{J}}_i$, respectively. Then, the model

$$\begin{array}{llll}\underset{\mathbf{x},\mathbf{y},\mathit{\lambda }}{\text{min}}& {\mathbf{c}}^{\mathrm{\top }}\mathbf{x}+{\mathbf{d}}^{\mathrm{\top }}\mathit{\lambda }& & \\ \text{s.t.} & \mathbf{y}\in {\text{Proj}}_{\mathbf{y}}\left(\text{NF}\right({\mathcal{D}}^{{\mathcal{J}}_i}\left)\right),& & \forall {\mathcal{J}}_i\in \{{\mathcal{J}}_1^1,\dots ,{\mathcal{J}}_k^1\}\\ & \mathbf{y}\in {\text{Proj}}_{\mathbf{y}}\left(\text{NF}\right({\mathcal{D}}_{\mathtt{\text{outer}}}^{{\mathcal{J}}_i}\left)\right),& & \forall {\mathcal{J}}_i\in \{{\mathcal{J}}_1^2,\dots ,{\mathcal{J}}_q^2\}\\ & {\mathit{T}}^{\mathrm{\top }}\mathit{\lambda }=\mathbf{y},\hspace{0.17em}\mathbf{y}\in \text{Relax}\left(\mathcal{Y}\right),\hspace{0.17em}\mathbf{y}\in \mathcal{S}\left(\mathbf{x}\right),\hspace{0.17em}\mathbf{x}\in \mathcal{X}& & \end{array}$$(11)

generates a feasible solution $\mathbf{x}\in \mathcal{X}$ and lower bound on the optimal objective value of model (1).

As a proof of concept, consider the model from Example 2 but with an additional constraint of $y_1+y_2+y_3\le 2$. Rather than formulating an exact or relaxed decision diagram of

we could generate two exact or relaxed decision diagrams, one for each constraint, and combine them into the multinetwork flow model (11).

In practice, multinetwork approximations can be useful when single-network representations are too large or when the underlying problem structure does not admit an obvious recursive formulation. In these settings, we can generate exact or relaxed DDs for subsets of constraints that do admit a simple recursive formulation (Tjandraatmadja and van Hoeve 2019, Davarnia and Van Hoeve 2021). Furthermore, in many decision-making problems, a large subset of constraints defining $\mathcal{Y}$ may have a totally unimodular constraint coefficient matrix and integer right-hand sides and thus define a feasible space for which its relaxation is an integral polytope. In this setting, we can generate DDs only for the remaining constraints. The corresponding multinetwork flow approximations may potentially be smaller in size, easier to implement, or better in quality than single-network approximations. We will further explore these models in the numerical experiments in Section 6.2.

5.3. Evaluating the Quality of a Solution

All models presented earlier generate a feasible first-stage solution $\widehat{\mathbf{x}}\in \mathcal{X}$. The true objective value of this solution, which we denote using $\text{OV}(\widehat{\mathbf{x}})$, can be computed by solving

Problem (12) can be solved using a simple constraint generation algorithm where, at an iteration k, a master problem solves a relaxation of (12) over an enumerated set of recourse solutions ${\widehat{\mathcal{Y}}}^k\subset \mathcal{Y}\cap \mathcal{S}(\mathbf{x})$, whereas a subproblem solves the inner problem for a fixed ${\widehat{\mathit{\xi }}}^k$ and generates a new solution ${\widehat{\mathbf{y}}}^k\in \mathcal{Y}\cap \mathcal{S}(\widehat{\mathbf{x}})$ to add to ${\widehat{\mathcal{Y}}}^k$; for brevity, we refer to Kämmerling and Kurtz (2020) for details.

Because the approximation models presented in Section 5.2 produce both a solution $\widehat{\mathbf{x}}\in \mathcal{X}$ and a dual bound, we can compute a model-based optimality gap for the solution, defined as

where “dual bound” denotes the optimal objective value of the approximation model. This model-based optimality gap provides an upper bound on the true optimality gap and allows us to quantify the quality of any solution $\widehat{\mathbf{x}}\in \mathcal{X}$. This is important both for practical purposes and for the numerical experiments, where the true optimality gap in large-scale ARBO problems cannot be computed because the underlying problems are too difficult to solve.

6. Numerical Experiments

In this section, we examine the performance of the exact and approximate network flow models for two different ARBO problems. All experiments and models were coded in Python 3.7 and solved using Gurobi 9.1.1 under default settings on a Macbook Pro (M1 Chip) with 16GB of RAM. The code can be found at Bodur et al. (2025).

6.1. Capital Budgeting Problems

We first consider a capital budgeting problem, which seeks to compute a robust investment plan for a set of n projects under an investment budget constraint. Several variants of this problem have been studied in the adaptive robust optimization literature (Hanasusanto et al. 2015, Kämmerling and Kurtz 2020, Subramanyam et al. 2020, Arslan and Detienne 2022, Dumouchelle et al. 2023). In this problem, investment decisions are binary and can be made in two stages; that is, projects can be invested in at either an early stage (denoted using $\mathbf{x}\in {\{0,1\}}^n$) or a later stage (denoted using $\mathbf{y}\in {\{0,1\}}^n$). The profitability of a project is not known with certainty in the first stage and is revealed only in the second stage. Early-stage investors are rewarded with a first-mover advantage that entails a higher percentage on the final profit generated by each project.

We follow the model formulation from Arslan and Detienne (2022) and use their publicly available problem instances.¹ Specifically, we consider the adaptive capital budgeting problem:

In this model, $g_i$ denotes the cost of investing in project i, h is the total investment budget, $f\in [0,1)$ captures the first-mover advantage, and ${\xi }_i$ is the payoff of project i, which is unknown in the first stage. The payoff of each project depends on a set of M common risk factors ${\alpha }_1,\dots ,{\alpha }_M$. Specifically, it is assumed that ${\xi }_i={\sum }_{j=1}^M{U}_{ij}{\alpha }_j$, where $U_{ij}\in \mathbb{R}$ describes the impact that each risk factor ${\alpha }_j$ has on the payoff ${\xi }_i$. The risk factors ${\alpha }_j$ are assumed to reside in $[-1,1]$, that is,

We consider 300 instances of varying size where $n=\{10,20,30,40,50\}$ and where each value of n is associated with 60 different instances (Arslan and Detienne 2022). Each instance is characterized by a random sample of project costs and payoffs as well as an investment budget defined as a fraction of total project costs, that is, $h=m{\sum }_{i=1}^n{g}_i$ with $m=\{0.2,0.4,0.6,0.8\}$.

In (14), the linking constraint $\mathbf{y}\ge \mathbf{x}$ satisfies selective adaptability, and the remaining two constraints can be reformulated using DDs. Because ${\mathbf{g}}^{\mathrm{\top }}\mathbf{y}\le h$ is a knapsack constraint, we generate the exact decision diagram by following the steps in Example 1; from this, we can obtain an exact network flow reformulation of (14). We also derive approximate network flow models based on relaxed DDs generated using the distance-based merging approach discussed in Section 5.1.2. Recall that in this approach, a user-specified value Q bounds the distance between state values of nodes that are to be merged. In the capital budgeting problem, the state value of a node in layer i defines the amount of weight that has been added to the knapsack based on decisions $y_1,\dots ,y_i$, and Q represents the maximum difference in weight values between any two nodes to be merged. We follow the procedure described in Example 4 to determine which nodes should be merged. Once merged, the new node is assigned the minimum state value of the nodes that were merged. Finally, recall that the resulting network flow models generate first-stage solutions and dual bounds.

6.1.1. Model Formulation and Solution Time.

Table 1 highlights the average number of arcs in the exact ($Q=0$) and approximate ($Q\ge 1$) decision diagrams used to represent ${\mathbf{g}}^{\mathrm{\top }}\mathbf{y}\le h$. These numbers convey the size of the corresponding network flow models because each arc in a decision diagram corresponds to a continuous variable in the model. All decision diagrams took no more than a few seconds to generate, and most were generated within a fraction of a second.

Table 1. The Average Number of Arcs in the Reduced Decision Diagrams Under Different Q Values

Table 2 summarizes the solution time of the network flow models, from which we make two observations. First, when n increases, the time it takes to solve the exact model grows exponentially, implying that the exact models are generally intractable beyond small values n. In contrast, as we increase Q for a fixed value of n, the solution time can be reduced drastically. For example, when $n=50$, the solution time of the models can be reduced by several orders of magnitude when $Q\ge 3$.

Table 2. The Average Solution Time (in Seconds) of the Network Flow Models Under Different Values of Q

6.1.2. Quality of Model Solutions.

For each solution $\widehat{\mathbf{x}}\in \mathcal{X}$ generated by the network flow model, we calculate its model-based optimality gap as well as the true optimality gap where possible (see Section 5.3). Table 3 summarizes these values. Note that only the model-based optimality gap is shown for the $n=50$ instances because the exact models (i.e., where $Q=0$) could not be solved within the time limit. We also note that because the solution times of the exact models are generally negligible for $n=10$ and $n=20$ (see Table 2), we focus our discussion on the instances with $n=30,40$ and 50, where tractable approximations become more critical.

Table 3. The Average True Optimality Gap and Model-Based Gap (Shown in Parentheses) of the Solutions

The main takeaway from Table 3 is that both the true and the model-based optimality gaps remain very small, even with the significant reduction in solution times shown in Table 2. For example, when $n=40$ and $Q=5$, the average optimality gap of solutions is less than 0.5% despite the models taking two orders of magnitude less time to solve compared with the exact model (i.e., 3 versus 534 seconds; Table 2). Similarly, when $n=50$, the average model-based optimality gap is $0.08\%$ when $Q=1$ and $0.5\%$ when $Q=10$, the latter of which requires an average solution time that is three orders of magnitude less than the exact model (i.e., 2 versus 3,600+ seconds; Table 2).

These observations highlight that the approximate network flow models are able to simultaneously generate both (i) near-optimal solutions and (ii) near-exact dual bounds. In Section EC.3.3 of the Electronic Companion, we provide a more in-depth analysis of the dual bounds and also compare our results to the K-adaptability approach from Hanasusanto et al. (2015) as well as the branch-and-price approach from Arslan and Detienne (2022).

6.2. Robust Assignment Problems

In this section, we consider a problem setting in which there are numerous linking constraints between first-stage and second-stage decisions. Specifically, we examine robust and adaptive assignment problems and focus on the use of multinetwork flow models (discussed in Section 5.2).

Assignment problems encompass numerous decision-making tasks that span many applications. A standard assignment problem can be modeled as a bipartite graph $(\mathcal{L},\mathcal{M},\mathcal{S})$ of agents $\mathcal{L}=\{1,\dots ,L\}$, tasks $\mathcal{M}=\{1,\dots ,M\}$, and directed links $\mathcal{S}$. Let $\mathcal{L}(m)\subseteq \mathcal{L}$ denote the subset of agents for which there exists a directed link to task $m\in \mathcal{M}$, and similarly, let $\mathcal{M}(\ell )\subseteq \mathcal{M}$ denote the subset of tasks that can be assigned to agent $\ell \in \mathcal{L}$. Depending on the context, agents can represent people, products, resources, funding, or jobs, whereas tasks can represent locations, facilities, projects, or machines. Each agent $\ell \in \mathcal{L}$ is associated with weight $a_{\ell }\ge 0$, and each task m has capacity $b_m\ge 0$, whereas the reward ${\xi }_{\ell ,m}\ge 0$ (e.g., match quality) of a specific agent-task pairing is unknown (e.g., varies day to day or must be estimated from data). Figure 2 illustrates this assignment problem.

Using this notation, recall that a static robust assignment problem can be formulated as

where $\mathcal{X}=\mathcal{Y}$ and $\mathcal{Y}$ is defined as

We now consider a setting where a planner desires a limited degree of flexibility in the assignment planning process by pre-identifying a subset of potential agent-task pairings from which assignments can be made once rewards become known. Practically, this may arise when there is a need to train people for specific tasks, inform individuals about potential assignments, or inspect match quality before making a final decision. We assume that we have a cardinality constraint limiting the number of preidentified pairings and formulate the adaptive robust assignment problem as

where $\mathcal{X}≔\{\mathbf{x}\in {\{0,1\}}^{|\mathcal{S}|}:{\Vert \mathbf{x}\Vert }_1\le \beta ·\left|\mathcal{S}\right|\}$. In our experiments, we consider $\beta$ values of 0.5, 0.6, 0.7, 0.8, and 0.9, which correspond to the ability to pre-identify $50\%$, $60\%$, $70\%$, $80\%$, and $90\%$ of pairings in the first stage, which can then be used to form a final assignment in the second stage.

6.2.1. Experimental Setup.

We consider assignment problems of five sizes where $(L,M)$ is equal to (20,2), (20,3), (20,4), (25,5), and (25,8). For each problem size, we randomly generate 10 instances. Each instance is characterized by (i) a randomly generated bipartite graph with 50% sparsity (such that $\left|\mathcal{S}\right|=20,30,40,63$ and 100, respectively) and (ii) randomly generated vectors of agent weights $\mathbf{a}$ and task capacities $\mathbf{b}$, where each element of $\mathbf{a}$ is a random integer from [1, 10] and each element of $\mathbf{b}$ is a random integer from $[2·\text{max}(\mathbf{a}),{\Vert \mathbf{a}\Vert }_1/m]$.

For each instance, we generate the uncertainty set as follows. Let ${\mathit{\xi }}^0$ denote a vector representing the nominal reward of each link. With slight abuse of notation, each element ${\xi }_i^0$ is a random number generated from [0.5 a, a], where a is the weight of the agent associated with this link. Then, the uncertainty set is defined on the percentage of deviation from ${\mathit{\xi }}^0$, namely,

We consider both an exact network flow formulation and a multinetwork flow approximation, abbreviated as Exact NF and Multi NF:

• Exact NF. Model (15) can be represented as an exact network flow model by reformulating $\mathcal{Y}$ based on a straightforward extension of the procedure outlined in Example 1. Specifically, the state $\mathit{S}$ in the recursive formulation of $\mathcal{Y}$ is a vector of size $L+M$ (rather than scalars) that captures the amount of capacity remaining in the constraints of $\mathcal{Y}$ based on previous decisions of $y_1,\dots ,y_i$. Details are given in Section EC.4.1 of the Electronic Companion.

• Multi NF. To form the multinetwork flow model, we generate an exact decision diagram for each knapsack constraint ${\sum }_{\ell \in \mathcal{L}(m)}{a}_{\ell }{y}_{\ell ,m}\le b_m$; we leave the constraints ${\sum }_{m\in \mathcal{M}(\ell )}{y}_{\ell ,m}\le 1$ as is because they satisfy the integral polyhedron property (see Proposition 4). The resulting model generates a feasible first-stage solution and a dual bound.

The Multi NF model defined above can be considered as one of the most natural or straightforward approximations of the ARBO problem, and we will show that it performs well in the numerical experiments. Nonetheless, we remark that there are many possible variations of this model that we could use to tighten or loosen the approximation. For example, instead of generating one exact decision diagram for each knapsack constraint, we could generate one for each pair of constraints, which would tighten the approximation quality but potentially lead to larger formulations. On the other hand, we could also generate an approximate decision diagram for each knapsack constraint (like in Section 6.1), which leads to a worse approximation but could significantly reduce the size and solution time of the model. In summary, various Multi NF models could be defined that will trade off between lower solution times and better solution quality.

Finally, for each network flow model, we enforce a one-hour time limit to build and reduce each diagram. Similar to Section 6.1, we also consider the K-adaptability model with $K=3$ and $K=4$ as a point of reference; the complete MILP formulation can be found in Section EC.4.3 of the Electronic Companion. For all models, we enforce a 30-minute time limit for the solver.

6.2.2. Computational Results.

Table 4 provides an overview of the size and the time taken to generate and reduce the DDs. The size and time for exact DDs increases rapidly with the size of the assignment problem; on average, it takes one second for the (20,2) instances, 73 seconds for the (20,3) instances, and more than one hour for all the other instances. Because the set of all solutions in $\mathcal{Y}$ must be represented by an exact DD, the size of this DD may grow exponentially with the number of constraints used to define $\mathcal{Y}$. In contrast, the average size and time for the multinetwork models is far smaller and generally negligible (<1 second). Because the Multi NF model is simply a collection of individual DDs that are each exact reformulations of only one constraint in $\mathcal{Y}$, the size and formulation time scales linearly in the number of constraints defining $\mathcal{Y}$.

Table 4. The Build Time and Attributes of the Decision Diagrams

Table 5 highlights the average solution time of the models over different problem instances. We observe that the Multi NF models can be solved relatively efficiently even as the size of the assignment problem grows larger. This is in direct contrast to the K-adaptability model.

Table 5. Average Solution Times Across the Different Models

Figure 3 shows the model-based optimality gap of solutions generated by the Multi NF models; as we discussed in Section 5.3, these values are upper bounds on the true optimality gap of the solution. First, we note that 86 out of the 250 problem instances had an optimality gap of 0%. In these instances, the Multi NF model found the optimal solution and verified its optimality. Most of these cases pertained to problem instances that were smaller in size. Second, 223 out of 250 instances had an optimality gap that was less than 1%, 247 out of 250 instances had a gap that was less than 2%, and no solution had an optimality gap that was greater than 3%. These statistics include solutions of the Multi NF instances that could not be solved within the time limit.

6.3. Summary of Numerical Results

We briefly summarize the main takeaways of the numerical section. First, we show that exact network flow formulations can be tractably formulated and solved for smaller instances of ARBO problems. Second, for larger instances, the approximate network flow formulations are tractable and simultaneously generate (i) near-optimal first-stage solutions and (ii) high-quality dual bounds. Finally, the size and solution times of the approximate network flow models can often be reduced drastically while sacrificing very little in terms of the quality of solutions and dual bounds.

7. Conclusion

In this paper, we examine ARBO problems with objective uncertainty. We leverage ideas from the decision diagram community to reformulate and approximate these problems as constrained network flow models, which can be solved using standard solvers. Our numerical experiments show that they can efficiently generate near-optimal solutions and high-quality bounds.

By forming this connection between robust optimization and decision diagrams, our framework can also take advantage of independent contributions that emerge from the decision diagram research community. For example, recent literature has explored the use of DDs to solve deterministic integer and/or nonlinear optimization problems (see Section 2). This suggests that it may be possible to extend our framework to more general adaptive robust optimization problems.