{'fraction_non_alphanumeric': 0.04743918551332867, 'fraction_numerical': 0.015357115100645061, 'mean_word_length': 4.123108775282688, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constrained Shortest-Path Reformulations for Discrete Bilevel and Robust Optimization

Many discrete optimization problems are amenable to constrained shortest-path reformulations in an extended network space, a technique that has been key in convexification, bound strengthening, and search mechanisms for a large array of nonlinear problems. In this paper, we expand this methodology and propose constrained shortest-path models for challenging discrete variants of bilevel and robust optimization, where traditional methods (e.g., dualize-and-combine) are not applicable compared to their continuous counterparts. Specifically, we propose a framework that models problems as decision diagrams and introduces side constraints either as linear inequalities in the underlying polyhedral representation, or as state variables in shortest-path calculations. We show that the first case generalizes a common class of combinatorial bilevel problems where the follower's decisions are separable with respect to the leader's decisions. In this setting, the side constraints are indicator functions associated with arc flows, and we leverage polyhedral structure to derive an alternative single-level reformulation of the bilevel problem. The case where side constraints are incorporated as node states, in turn, generalizes classical robust optimization. For this scenario, we leverage network structure to derive an iterative augmenting state-space strategy akin to an L-shaped method. We evaluate both strategies on a bilevel competitive project selection problem and the robust traveling salesperson with time windows, observing considerable improvements in computational efficiency as compared to current state-of-the-art methods in the respective areas.


Introduction
Given a graph G equipped with arc lengths, the constrained shortest-path problem (CSP) asks for the shortest path between two vertices of G that satisfies one or more side constraints, such as general arc-based budget restrictions or limits on the number of nodes traversed. CSPs are fundamental models in the operations literature, with direct applications both in practical transportation problems (Festa 2015), as well as in a large array of solution methodologies for routing and multiobjective optimization (Irnich and Desaulniers 2005). Due to its importance, the study of scalable algorithms for CSPs and related problems is an active and large area of research (e.g., Cabrera et al. 2020, Vera et al. 2021, Kergosien et al. 2022.
In this paper, we expand upon this concept and propose the use of CSPs as a modeling construct for more general classes of optimization problems. Our approach consists of reformulating either a relaxation or a projection of the problem in an extended compact network space, where paths of the network have a one-to-one relationship to potential feasible solutions of the system. Solving such a relaxation or projection then reduces to finding a constrained shortest path, where side constraints incorporate information of the missing constraints or variables of the original problem.
In particular, we exploit network reformulations based on decision diagrams (Bergman et al. 2016), which are compressed networks obtained from the state-transition graph of recursive models.
The distinct attribute of this approach is that one can leverage network structure to derive new reformulations and more efficient algorithm classes, which provide a novel "combinatorial" perspective to such problems when the framework is applicable. More precisely, we exploit the dual description of a constrained shortest-path problem as either a network-flow model (i.e., its polyhedral perspective), or as a label-setting search (i.e., its dynamic programming perspective), applying each based on the setting at hand and the information captured by the side constraints.
Specifically, in this work we focus and illustrate the modeling construct on two classes of optimization problems. The first refer to binary bilevel programs of the form where variables x F and x L formulate a leader's and a follower's decisions, respectively. The model above is prevalent in network interdiction problems (Morton et al. 2007, Cappanera and Scaparra 2011, Hemmati et al. 2014, Lozano and Smith 2017a, minimum edge and vertex blocker problems (Bazgan et al. 2011, Pajouh et al. 2014, and other classes of adversarial settings (Costa et al. 2011, Caprara et al. 2016, Zare et al. 2018). However, its challenge stems from the possibly difficult combinatorial structure associated with the follower's variables, often preventing direct extensions of duality-based solution methodologies from continuous bilevel programs.
We show that, alternatively, the follower's subproblem can be formulated as a CSP, with univariate side constraints parameterized by the leader's variables x L . In particular, the coefficient matrix of the CSP has a totally unimodular structure and, thus, allows typical polyhedral approaches (i.e., dualize-and-combine) to rewrite the "argmax" of the follower's problem as a feasibility system. Further, we also exploit the property of the duals of the follower's network reformulation to convexify non-linear inequalities associated with this model, which appear due to complementary slackness constraints derived from the dualize-and-combine technique. The resulting model is an extended, single-level mixed-integer linear programming (MILP) reformulation of the original bilevel problem, and thus amenable to commercial state-of-the-art solvers.
The second problem class we investigate is robust optimization, min x∈{0,1} {cx : where ∆ is a uncertainty set that parameterizes the realizations of the coefficient matrix A(δ) and right-hand side b(δ). Problems of this form have been investigated extensively, with existing approaches primarily based on cutting-plane algorithms akin to Benders decomposition (Mutapcic and Boyd 2009, Zeng and Zhao 2013, Ben-Tal et al. 2015, Ho-Nguyen and Kılınç-Karzan 2018, Borrero and Lozano 2021. Each step in such procedures involves solving an MILP of increasing size, which may often inhbit the scalability of the approach. The second CSP model we consider is based on dynamic programming (DP) and related models, which have been used extensively in state-of-the-art techniques for large-scale CSPs (e.g., Dumitrescu and Boland 2003, Cabrera et al. 2020, Vera et al. 2021. In general, the DP in this context solves a shortest-path problem on an extended graph via a specialized recursion that labels the states achievable at a node. The benefit of such models is that they do not require linearity and exploit the combinatorial structure of the graph for efficiency. We expand on such models in §4 and refer to Irnich and Desaulniers (2005) for general survey of labeling algorithms for CSPs.
Our first main contribution is a reformulation of a discrete class of bilevel optimization problems. Bilevel optimization is an active research field with applications in energy and natural gas market regulation (Dempe et al. 2011, Kalashnikov et al. 2010, waste management (Xu and Wei 2012), bioengineering (Burgard et al. 2003), and traffic systems (Brotcorne et al. 2001, Dempe andZemkoho 2012), to name a few. Variants where the follower's problem is convex are solved by a dualize-and-combine approach, which combines the Karush-Kuhn-Tucker (KKT) optimality conditions of the follower's problem with the leader's model. This results in a single-level nonlinear reformulation approachable by a wide array of convexification and other specialized approaches; we refer to the survey by Dempe et al. (2015), Kleinert et al. (2021) for further details and applications.
The dualize-and-combine approach, however, is not generally applicable for cases where decision variables are discrete. An important special case refers to the large class of interdiction problems, where the leader and follower play an adversarial relationship (see, e.g., survey by Smith and Song 2020). Discrete bilevel problems are generally challenging and have fostered a large array of methodologies that extended their single-level counterpart. Solution approaches include branch and bound (DeNegre andRalphs 2009, Xu andWang 2014), Benders decomposition (Saharidis and Ierapetritou 2009), parametric programming (Domínguez and Pistikopoulos 2010), column-and-rown generation (Baggio et al. 2021), and cutting-plane approaches based on an optimal value-function reformulations (Mitsos 2010, Lozano andSmith 2017b). In particular, the state-of-the-art approach for mixed-integer linear bilevel programs is the branch-and-cut algorithm by Fischetti et al. (2017), which we use as a benchmarck algorithm in our computational experiments. In contrast to such literature, our method is a dualize-and-combine technique that is applicable when the leader wishes to "block" the follower's decisions. We exploit the structure of the decision diagram to provide a linear optimality certificate for the follower's problem, which leads to an alternate single-level reformulation that is amenable to classical mathematical programming solvers.
Our second contribution is an approach to address discrete robust optimization problems via CSP reformulations. Robust optimization is also pervasive in optimization, with extensive literature in both theory , Ben-Tal et al. 2006, Bertsimas and Brown 2009, Li et al. 2011, Bertsimas et al. 2016) and applications (Lin et al. 2004, Ben-Tal et al. 2005, Bertsimas and Thiele 2006, Yao et al. 2009, Ben-Tal et al. 2011, Gregory et al. 2011, Moon and Yao 2011, Gorissen et al. 2015, Xiong et al. 2017. For settings where variables are discrete, state-of-the-art algorithms first relax the model by considering only a subset of realizations, iteratively adding violated variables and constraints from missing realizations until convergence. This primarily includes cutting plane algorithms based on Benders decomposition (Mutapcic and Boyd 2009, Zeng and Zhao 2013, Ben-Tal et al. 2015, Ho-Nguyen and Kılınç-Karzan 2018, Borrero and Lozano 2021. Such methodologies exploit, e.g., the structure of the uncertainty set to quickly identify violated constraints via a "pessimization" oracle that finds the worst-case realization of the uncertainty set. Our methodology can also be perceived as a type of cutting-plane approach that augments an uncertainty set initially composed of a small number of realizations. In particular, the approach incorporates the scenario-specific constraints as resources in a full CSP reformulation of the robust problem based on decision diagrams, which we later remodel as a infinite-dimensional dynamic program. We address such a program by adding missing state iteratively, where the new state captures one or more violated constraints of the original model. Existing state-augmenting algorithms have been applied by Boland et al. (2006) in labeling methods to solve the standard CSP, as well as a in dynamic programs for stochastic inventory management (Rossi et al. 2011).

Preliminaries
We introduce in this section the concept of decision-diagram (DD) reformulations that we leverage throughout this work. We start in §3.1 by introducing the network structure and notation. Next, we briefly discuss in §3.2 existing construction strategies and problems amenable to this model.

Network Reformulations via Decision Diagrams
In this context, a decision diagram D is a network reformulation of the problem where X ⊆ Z n is an n-dimensional finite set for n ≥ 1. Specifically, D = (N , A, ν, l) is an orderedacyclic layered digraph with node set N and arc set A ⊆ N × N . The set N is partitioned by n + 1 layers N 1 , . . . , N n+1 , where N 1 = {r} for a root node r and N n+1 = {t} for a terminal node t.
Each arc a ∈ A is equipped with a value assignment ν a ∈ Z and a length l a ∈ R. Further, arcs only connect nodes in adjacent layers, i.e., with each a ∈ A we associate a tail node t a ∈ N j and a head node h a ∈ N j+1 , j ∈ {1, . . . , n}. We denote by τ a ∈ {1, . . . , n} the layer that includes the tail node of a, i.e., t a ∈ N τa .
The decision diagram D models DO through its r − t paths and path lengths, as follows. Given an arc-specified r − t path (a 1 , . . . , a n ), where h a i = t a j+1 for j = 1, . . . , n − 1, the solution x = (ν a 1 , . . . , ν an ) composed of the ordered arc-value assignments is such that x ∈ X and the path length satisfies n j=1 l a j = f (x). Conversely, every solution x ∈ X maps to a corresponding r − t where its length matches the objective evaluation of x.
Thus, by construction, the shortest path in D yields an optimal solution to DO. We also remark that if DO were a maximization problem, the longest r − t path provides instead the optimal solution, which is also computable efficiently in the size of D since the network is acyclic.
Example 1. Figure 1(a) depicts a reduced decision diagram for the knapsack problem Castro et al. 2021). Since variables are binaries, the value assignment of an arc a ∈ A is either ν a = 0 (dashed lines) or ν a = 1 (solid lines). An arc a emanating from the j-th layer, i.e., t a ∈ N j , corresponds to an assignment x j = ν a . Further, the length of a is either 0 if ν a = 0 or the coefficient of x j in the objective otherwise. In particular, the longest path (in blue) provides the optimal solution x = (0, 0, 1, 1) with objective 15. Example 2. Green-shaded arcs represent the longest-path and optimal solution (in colour.)

Modeling Problems as DDs
The standard methodology to model DO as D consists of two steps (Bergman et al. 2016). The first reformulates DO as a dynamic program (DP) and generates its state-transition graph, where states and actions are mapped to nodes and arcs, respectively. The second steps compresses the graph through a process known as reduction, where nodes with redundant information are merged. This is a key process in the methodology, often reducing the graph by orders of magnitude in practice (e.g., Newton and Verna 2019). We detail the modeling framework in Appendix F for reference.

Constrained Shortest-Path Models
In this section, we introduce the two constrained shortest-path models over DDs that we leverage in our CSP reformulations. Specifically, we present a mixed-integer linear programming model (MILP) in §4.2 and a dynamic program (DP) in §4.3. The first model exposes polyhedral structure and is appropriate, e.g., if the resulting formulation is sufficiently compact for mathematical programming solvers. The second model preserves the decision diagram structure and allows for combinatorial algorithms exploiting the network for scalable solution approaches.

General Formulation
Let D = (N , A, ν, l) be a decision diagram encoding a given discrete optimization problem DO. We investigate reformulations of optimization problems based on the constrained shortest-path model min y∈{0,1} |A| a∈A l a y a (CSP-D) where the inequality system (2) corresponds to m ≥ 1 arc-based side constraints for a non-negative coefficient matrix G ∈ R m×|A| and a non-negative right-hand side vector d ∈ R m .
We consider additional modeling assumptions when reformulating problems as CSP-D. First, D is of computationally tractable size in the number of variables n of the model DO it encodes (e.g., as in the cases in §3.2). Second, building a new decision diagram D that also satisfies a subset of constraints (2) is not practical due to, e.g., its potentially large size. Finally, for simplicity we also assume that there exists at least one r − t path that satisfies (2). Results below can be directly adapted otherwise to detect infeasibility.
Our goal is to study CSP-D as a modeling approach to reveal structural properties and derive alternative solution methods to broader problem classes. In particular, modeling problems as CSP-D consists of identifying which constraints should be incorporated within D, and which should be cast in the form of side constraints (2) by choosing an appropriate G and d.
We note that CSP-D is a special case of the resource-constrained shortest-path problem (Irnich and Desaulniers 2005) on the structure imposed by a decision diagram. For instance, if variables are binaries, each node has at most two outgoing arcs. Nonetheless, model CSP-D still preserves generality and hence it is of difficult computational complexity, as we indicate below.
Proposition 1. CSP-D is (weakly) NP-hard even when m = 1 and D has one node per layer.

MILP Model
is a binary linear reformulation of CSP-D which is suitable, e.g., to standard MILP solvers.
The objective function of MILP-D, alongside constraints (3)-(4), is a perfect linear extended formulation of DO in the space of y variables (Behle 2007). That is, there exists a one-to-one mapping between solutions y satisfying (3)-(4) and a feasible solution x ∈ X , and the objective function value of y match that of f (x). Thus, model MILP-D provides a polyhedral representation of the constrained shortest-path problem even when X contains non-trivial combinatorial structure, which is captured by D and that will be exploited in our reformulations. We also note that similar models combining a shortest-path linear formulation with side constraints appear, e.g., in classical extended formulations based on DPs (e.g., Eppen andMartin 1987, Conforti et al. 2010).

Dynamic Programming Model
An alternative common model is to reformulate CSP-D as a DP, which allows for combinatorial algorithms that exploit the network structure and scale to larger problem sizes. The underlying principle is to perceive each inequality in (2) as a "resource" with limited capacity. DP approaches enumerate paths ensuring that the consumed resource, here modeled as state variables, does not exceed the budget specified by the inequality right-hand sides.
We introduce such a DP model for CSP-D defined over the structure of a decision diagram D. With each node u ∈ N we associate a function V u : R m → R that evaluates to the constrained shortest-path value from u to the terminal t. The function argument is a state vector s = (s 1 , . . . , s m ) ∈ R m where each s i represents the amount consumed of the i-th resource, i ∈ {1, . . . , m}.
Let G a be the vector defined by the a-th column of G, a ∈ A. We write the Bellman equations where min ∅ {·} = +∞. Note that DP-D resembles a traditional shortest-path forward recursion in a graph (Cormen et al. 2009), except that traversals are constrained by the condition s + G a ≤ d.
Proof of Proposition 2. We will show this result by presenting a bijection between solutions of CSP-D and V r (0), also demonstrating that their solution values match. Specifically, for each node u ∈ N , let π u : R m → A be a function that maps a state s ∈ R m to any arc a ∈ A that minimizes V r (s), i.e., a policy. Since a ∈ Γ + (u) in the minimizer of V u (·) for any node u, unrolling the recursion V r (0) in DP-D results in a r − t path (u 1 , u 2 , u 3 , . . . , u n+1 ) such that u 1 = r, u n+1 = t, and whereŝ 1 = 0 andŝ j =ŝ j−1 + G πu j−1 (ŝ j−1 ) for j = 2, . . . , n. Thus, the evaluation of V r (0) matches the length of the path (u 1 , u 2 , u 3 , . . . , u n+1 ) considering the arcs specified by the policy π. Moreover, because of the condition s + G a ≤ d in the minimizer of V u (·) for any u ∈ N , we must haveŝ j ≤ d for all j = 1, . . . , n; in particular, for j = n we obtain for y a = 1 if π u j (ŝ j−1 ) = a for some j ∈ {1, . . . , n}, and 0 otherwise. That is, every path in D is feasible to CSP-D. Conversely, every solution to CPSP can be analogously represented as a policy π constructed with the arcs of the associated r − t path, noting that G non-negative implies that every partial sum of arcs also satisfy s + G a ≤ d.
Insights into the distinction between MILP-D and DP-D are revealed by duality. Given Proposition 2, let S u ⊆ R m be the set of states that are reachable at a node u ∈ N by the recursion starting with V r (0) and budget d, i.e., S r = {0} and The set S u is finite for all nodes u because D has a finite number of arcs. We can therefore reformulate DP-D as a linear program with the evaluations of V u (·) as variables: Notice that, because of (8), inequalities (7) simplify after replacing all variables V ha (·) by 0 whenever h a = t. The dual of the linear program obtained after such an adjustment is Thus, DP-D solves a (non-constrained) shortest-path problem on an extended variable space w a,s that incorporates the state information s, as opposed to the DD arc-only based variables y a in MILP-D. A solution to DP-D is equivalently a state-arc trajectory (s 1 , a 1 ), (s 2 , a 2 ), . . . , (s n , a n ) such that (a 1 , . . . , a n ) is an r − t path in D and s i = s i−1 + G a i−1 ≤ d for i = 2, . . . , n, s 1 = 0.
Example 2. Consider the knapsack instance from Example 1 and the side constraint We can formulate it as CSP-D using the DD from Figure 1(a) to encode the original feasible set X , plus a single side constraint a∈A g a y a ≤ 7 where g a is the coefficient of variable x τa in (11) if ν a = 1, and 0 otherwise. Figure 1(b) depicts the extended network space where LP-D is defined.
The numbers on the node circles represent the reachable states (i.e., labels) and the squares the original nodes in Figure 1(a) they are associated with. There exists a single label since m = 1.
Further, two nodes are connected if the labels are consistent with the side constraint and the arc exists in the original DD. For example, the labels in u 4 only have outgoing arcs with value 0, and their label remain constant since the coefficient of such arcs is also zero.
State-of-the-art combinatorial methods for general constrained path problems, such as the pulse method (Cabrera et al. 2020) and the CHD (Vera et al. 2021), are directly applicable and operate on the network specified in LP-D by carefully choosing which states to enumerate, a process commonly known as labeling. Such methods are effective if such labels have exploitable structure, which we leverage in §6. Otherwise, if the set of states is difficult to describe but the representation is sufficiently compact, then MILP-D could be more beneficial, which we investigate in §5.

Discrete Bilevel Programs
In this section, we introduce single-level reformulations based on CSP-D for bilevel programs where hand side vectors of the leader's and follower's, respectively; and 1 is an n-vector of ones. We assume the relatively complete recourse property, i.e., for any feasible leader solution there exists a corresponding feasible follower response, ensuring the existence of an optimal solution for the follower problem described in (13).
Formulation DB is prevalent in combinatorial bilevel applications from the literature. Specifically, the inequality x F ≤ 1 − x L in (13) represents exogenous "blocking" decisions by a leader to prevent actions from an non-cooperative follower, which operates optimally according to its own utility function. Examples of problems modeled as DB include network interdiction problems (Morton et al. 2007, Cappanera and Scaparra 2011, Hemmati et al. 2014, Lozano and Smith 2017a, minimum edge and vertex blocker problems (Bazgan et al. 2011, Pajouh et al. 2014, and other classes of bilevel adversarial problems (Costa et al. 2011, Caprara et al. 2016, Zare et al. 2018).
Problems of this class, however, are notoriously difficult because of the "argmax" constraint (13) that embeds the follower's optimization problem into the leader's model. In particular, the follower problem is non-convex since x F are binaries, preventing the use of typical techniques for continuous problems such as dualize-and-combine (Kleinert et al. 2021). Existing state-of-the-art methodologies rely, e.g., on specialized cutting-plane and search techniques (Fischetti et al. 2017).
In this section, we present a reformulation of DB as MILP-D, which is suitable to standard mathematical programming solvers if the underlying network D is of tractable size. We start in §5.1 by rewriting constraint (13) as a constrained shortest-path model, and extract an optimality certificate via polyhedral structure. Next, we present in §5.2 the full MILP reformulation based on the polyhedral structure of MILP-D. Finally, we discuss in §5.4 generalizations of the technique for broader classes of bilevel programs.

Reformulation of the Follower's Optimization Problem
Let D F = (N , A, ν, l) encode the discrete optimization problem DO with feasible set and objective f (x) = n j=1 c F j x j . That is, paths in D F correspond to solutions that satisfy the follower's constraints A F x ≤ b F and have arc lenghts given by Proposition 3 shows the link between constraint (13) and CSP-D.
Proposition 3. Let x L ∈ {0, 1} n be any feasible leader's solution in DB. There exists a one-to-one mapping between a follower's solution x F satisfying (13) and an optimal solution of CSP-D defined over D F with a maximization objective sense and m ≤ |A| side constraints of the form That is, the columns of G associated with arcs a ∈ A such that ν a = 1 form an identity matrix, and the columns for the remaining arcs are defined by zero vectors. Further, d a = 1 − x L τa if ν a = 1 and d a = 0 otherwise.
Proposition 3 leads to a combinatorial description of solutions x F observing (13) as paths in D F subject to side constraints (16). Thus, we can derive an equivalent mathematical program reformulation of constraint (13) through MILP-D. Inequalities (16), however, have a special form in that they only impose bounds on arc variables y, not significantly changing the coefficient matrix structure. Based on this observation, we present a more general result applicable when side constraints (2) do not break the integrality of the network-based extended model (3)-(4) presented in §4.2. We recall that G a is the a-th row of G, a ∈ A.
Proposition 4. Let D = (N , A, ν, l) encode a discrete optimization problem DO and let N be the coefficient matrix of y associated with the linear system (3)-(4). If the matrix N G is totally unimodular and d is integral, then the set of optimal solutions to MILP-D is Proj y Y, where Proof of Proposition 4. The assumptions of the statement imply that the optimal basic feasible solutions of the linear program max y a∈A are optimal to MILP-D and vice-versa, where g i,a is the element at the i-th row and a-th column of G. Consider the dual obtained by associating variables π with constraints (22)-(23) and variables γ with (24), i.e., Then, by strong duality, any optimal solution to both systems satisfy the system Y composed by the primal and dual feasibility constraints, in addition to inequality (19) where π and γ are of appropriate dimensions in terms of Y. We note that, to reduce the number of binary variables, we enforced integrality on x F as opposed to y, which does not impact the feasibility set with respect to x F because of constraint (35). Further, the system (29)-(36) follows the structure of KKT-based reformulation approaches, i.e., (29)-(31) ensure primal feasibility, (32)-(33) require dual feasibility, and (34) ensures strong duality.

Leader's MILP Model
An extended, single-level reformulation of DB can be obtained directly when replacing (13) by (29)-(36). However, the resulting formulation is nonlinear because of the terms γ a (1 − x L τa ) in (34). We show next how to linearize (34) based on Proposition 5.
Proposition 5. Any feasible solution (x F ,ŷ,π,γ) to (29)-(36) can be modified to satisfyγ a = M a for all arcs a ∈ A such that x L τa = 1 and ν a = 1 preserving feasibility, where M a values are large constants that satisfy M a ≥ l a +π ha −π ta , ∀a ∈ A : x L τa = 1 and ν a = 1.
Proof of Proposition 5. Consider any arbitrary feasible solution (x F ,ŷ,π,γ) to (29)-(36). Pick any arc a ∈ A such that x L τa = 1, ν a = 1, butγ a = M a . Notice thatγ a is in constraints (32) and (34). For (34), the difference (1 − x L τa ) is zero, and hence changingγ a does not impact such an equality. For (32) we restrict our attention to the caseγ a > M a ; notice that, otherwise, increasingγ a would not affect the inequalities. In this case, condition (37) guarantees that (32) is satisfied if we set γ a = M a , completing the proof.
Using the result above, we rewrite the left-hand side of (34) as where the second equality follows from Proposition 5. This yields the following MILP reformulation of DB: The inequality (41) ensures consistency of γ in line with Proposition 5 and the derivation of (40). We remark that MILP-DB is a standard MILP where only variables x L and x F are binaries.
Thus, the formulation is appropriate to standard off-the-shelf solvers whenever D F is a suitable encoding of the follower's suproblem.
The model above contains big-M constants, which is typical in specialized bilevel techniques (Wood 1993, Cormican et al. 1998, Brown et al. 2005, Morton et al. 2007, Bayrak and Bailey 2008. We show in Proposition 6 that under some mild conditions, a valid value for M a is given by l a for all a ∈ A : ν a = 1 and present in Proposition 7 a general way of obtaining valid M -values for MILP-DB.
Proposition 6. If the follower constraint coefficient matrix A F ∈ R m F ×n + , then setting M a = l a for all a ∈ A : ν a = 1 yields a valid formulation.
Proof of Proposition 6 First consider that if c F j < 0 then follower variable x F j = 0 at any optimal solution of DB since A F ∈ R m F ×n + . Thus we assume without loss of generality that c F j ≥ 0 for j ∈ {1, . . . , n}, which ensures that l a ≥ 0 for all a ∈ A.
We now show that there exists an optimal solution to MILP-DB that satisfies γ a = l a for all arcs a ∈ A such that x L τa = 1 and ν a = 1. Consider an optimal solution to MILP-DB given by (x L * , x F * , y * , π * , γ * ) and a decision diagramD = (N ,Â, ν, l) that is obtained by removing from A any arcs for which x L * τa = 1 and ν a = 1, i.e., removing from D F all the arcs with zero capacity in constraint (31). Formally, the modified set of arcs is given bŷ It follows from the non-negativity of A F that every node in D F has at least one outgoing zero arc and thus the set of nodes ofD is the same as the original set of nodes.
Since we only remove arcs with zero capacity, any feasible path inD is also a feasible path in D F , under x L * . Thus, the optimal path given by y * over D F is also an optimal primal solution over D. Next, let us consider an optimal dual solutionπ, corresponding to the optimal path given by y * overD. Given that all the arcs inÂ have a capacity of 1, we remove constraints (31) from the primal, which in turn removes variables γ from the dual. Dual feasibility ofπ hence ensures that π ta −π ha ≥ l a ∀a ∈Â and strong duality ensures that l y * −π r = 0.
We now show thatπ ta −π ha ≥ 0 for all a ∈ A. Note that dual variablesπ capture the node potential and are computed as the length of an optimal completion path from any node in N to t. These completion paths exists for every node in N since A F is non-negative.
Pick any arc a ∈ A such that ν a = 0. Since a ∈Â and l a = 0, constraints (43) ensure that π ta −π ha ≥ 0. Now pick any arc a ∈ A such that ν a = 1. Assume by contradiction thatπ ta −π ha < 0.
Since A F is non-negative, there exists an arc a ∈Â leaving t a with ν a = 0. Furthermore, observe that any partial solution composed of paths starting from h a and ending at t are also found in some other path starting at h a and ending at t, again because A F is non-negative and fixing x F τ a = ν a = 0 can only increase the number of solutions. As a result we have thatπ h a ≥π ha . Since a ∈Â and l a = 0, constraints (43) ensure thatπ t a ≥π h a . Because t a = t a we obtain that π ta ≥π h a ≥π ha , which contradicts our assumption thatπ ta −π ha < 0.
Using the results above now consider for any optimal solution to MILP-DB given by (x L * , x F * , y * , π * , γ * ), an alternative solution given by (x L * , x F * , y * ,π, γ ), where γ a = l a if x L * τa = 1 and ν a = 1 0 otherwise ∀a ∈ A.

Case Study: Competitive Project Selection
We present a case study in competitive project selection, a general class of bilevel knapsack problems that includes, e.g., Stackelberg games with interdiction constraints (Caprara et al. 2016). In this setting, two competing firms select projects to execute from a shared pool in a sequential manner, starting with the leader. Both firms seek to maximize their own profit, which is a function of both the projects as well as the competitor's selection. Applications of this setting are pervasive in Marketing, where a company must choose advertisement actions that are also impacted by the actions of their competitors (DeNegre 2011).
Let P = {1, . . . , n} be a set of n projects in a pool shared by two firms. The leader firm considers both the net profit of projects executed, given by c L j , and a penalty for projects executed by the follower, denoted by d L j . The follower firm considers only the net profit of their executed projects, given by c F j . Similarly, the project costs and the firm budgets are a L j , a F j and b L , b F , respectively. We write the competitive project selection problem as the bilevel discrete program max x F ∈ arg max x L ∈ {0, 1} n .
The objective function of (CPSP) maximizes the leader's profit penalized by the follower's actions.
The knapsack inequality (46) enforces the leader's budget. The follower's subproblem, in turn, is represented by (47) and maximizes profit subject to an analogous budget constraint. Moreover, the constraint x F ≤ 1 − x L imposes that only projects not selected by the leader can be picked.

Single-Level Reformulation.
We construct D F = (N , A, ν, l) based on a standard DP formulation for the knapsack problem for which the state variable stores the amount of budget used after deciding if project j ∈ P is selected or not (see e.g., Bergman et al. (2016)). We generate D F from the state transition graph of the DP described above by removing nodes corresponding to the infeasible state and merging all terminal states.
Example 3. Figure 2 shows a decision diagram generated for a follower's problem with |P| = 3 having a F 1 = a F 2 = 2, a F 3 = 4, and b F = 5. Solid lines represent yes-arcs, dashed lines represent no-arcs, and contributions to the follower's profit are shown alongside the yes-arcs.

Figure 2
A decision diagram corresponding to a follower's problem with 3 projects.
Each layer of the diagram corresponds to a follower's variable. There are 5 paths in the diagram that encode all the follower selections of projects that satisfy the budget constraint.
Following the approach described in Section 5.2, we reformulate CPSP over D F as the following single-level MILP: The objective function is the same as in the original formulation. Constraint (49) corresponds to the leader's budget constraint. Constraints (50) remain exactly as defined in Section 5.2, ensuring that y corresponds to a feasible path in D F (primal solution), that (π, γ) is a feasible dual solution, enforcing the relationship between flow-variables y and follower variables x F , and ensuring that y and (π, γ) satisfy strong duality. Additionally, since all the project costs are nonnegative, we set M a = c F τa according to Proposition 6.

Numerical Analysis.
We compare CPSP-D, dubbed DDR for decision-diagram reformulation, against an state-of-the-art approach for linear bilevel problems by Fischetti et al. (2017) based on branch and cut, denoted here by B&C. Our experiments use the original B&C code provided by the authors. For consistency, DDR was implemented using the same solver as B&C (ILOG CPLEX 12.7.1). All runs consider a time limit of one hour (3,600 seconds). We note that, while B&C is based on multiple classes of valid inequalities and separation procedures, the DDR is a stand-alone model. We focus on a graphical description of the results in this sections; detailed tables and additional results are included in Appendix I. Source code is available upon request.
We generated instances of CPSP based on two parameters: the number of items n, and the right-hand side tightness t = , smaller values of t correspond to instances with fewer feasible solutions. We considered five random instances for each n ∈ {30, 40, 50} and t ∈ {0.10, 0.15, 0.20, 0.25}, generating 60 problems in total. Instances for each parameter configuration are generated as follows. The coefficients a L are drawn uniformly at random from a discrete uniform distribution U(1, 25) and then we set a F = a L . Budgets b L = b F are set according to t. The project profits are proportional to their costs by setting c L i = 5a L i + ξ L i and c F i = 5a F i + ξ F i , where ξ-values are drawn independently and uniformly at random from a discrete uniform distribution U(1, 10).   seconds on average by DDR. This is at least one order of magnitude faster than B&C, which could not solve the majority of instances with t ≥ 0.2 within the time limit. We note that an analogous analysis follows for Figure 3b, as larger values of n increase the difficulty of the problem for both methods and the same scaling benefits are observed.
Based on these results, we generate two new classes of challenging instances to assess at which point DDR would still be sufficiently compact to provide benefits with respect to B&C. To this end, the first and second class draw the knapsack coefficients uniformly at random from the discrete uniform distributions U(1, 50) and U(1, 100), respectively, as opposed to U(1, 25). Larger domain distributions significantly increase the DD sizes, since the number of nodes per layers is bounded above by such numbers (we refer to Appendix I for the final DD sizes). We also vary the tightness within the set {0.1, . . . , 0.9} which also impacts the total number of nodes as discussed above. The number of variables is fixed to n = 30, generating in total 45 instances for each new class with five random instances per configuration. (Coloured) Performance profile plot for varying uniform discrete distribution of the knapsack constraint coefficients a F and a L .

Generalizations
We end this section by discussing extensions of our approach to more general problem settings.
Consider a problem setting given by and objective f (x) = g F (x), i.e., paths in D F correspond to follower solutions that satisfy H F (x) ≤ 0 and path lengths correspond to exact evaluations of function g F . Our reformulation procedure and all the results in propositions 3 -7 can be directly extended for the more general setting above by defining the side constraints as because there is nothing in our methodology that requires the leader or follower problem to be linear. The main two conditions for our results to hold are that the follower problem is correctly encoded by D F and that the effect of the leader variables on the follower problem happens exclusively via linking constraints of the form x F i ≤ 1 − µ i (x L ). Using our approach, we generate a single-level reformulation of (51) given by x L ∈ {0, 1} n .
We remark that (54) is a single-level problem that conserves the original functions g L , H L , and µ i while replacing functions g F and H F with a series of linear constraints based on the variables of a primal flow problem over D F and the corresponding dual variables.

Discrete Robust Optimization
In this section, we introduce a CSP-D reformulation to robust problems of the form where X is a discrete set as in DO, δ are realizations of a (possibly infinite) uncertainty set ∆, and A(δ) ∈ R r×n , b(δ) ∈ R r are the coefficient matrix and the right-hand side, respectively, for the realization δ ∈ ∆ for some r > 0. As discussed in §2, state-of-the-art algorithms for DO first relax the model by considering only a subset of ∆ in (55), iteratively adding violated constraints from missing realizations until convergence.
We propose an alternative algorithm for RO that is applicable when the set X is amenable to a DD encoding. The potential benefit of the methodology is that it is combinatorial and exploits the network structure of D for scalability purposes. Our required assumptions are the existence of a separation oracle S (x) that identifies violations of (55), i.e., and that such oracle returns either a pair (i, δ) or ∅ in finite computational time. This is a mild assumption in existing RO formulations and holds in typical uncertainty sets, such as when ∆ is an interval domain or, more generally, has a polyhedral description. For instance, if ∆ is finite (e.g., derived from a sampling process), then the simplest S (x) enumerates each δ separately.
We start in §6.1 with a reformulation of the robust problem as a CSP-D, leveraging this time the dynamic programming perspective presented in §4.3. Next, we discuss a state-augmenting procedure to address the resulting model in §6.2, and perform a numerical study on a robust variant of the traveling salesperson problem with time windows (TSPTW) in §6.3.

Reformulation of RO
Let D R = (N , A, ν, l) encode the discrete optimization problem DO with feasible set X and objective f (x). The constrained longest-path reformulation is obtained directly by representing (55) as side constraints over the arc space of D R , which we formalize in Proposition 8.
Proposition 8. For each arc a of D R , δ ∈ ∆, and i ∈ {1, . . . , r}, define the scalars where a i,τa (δ) is the (i, τ a )-th element of A(δ). Let G(δ) = {g i,a (δ)} ∀i,a be the matrix composed of such scalars for each δ. There exists a one-to-one mapping between solutions x of RO and paths of CSP-D defined over D R and with side constraints Modeling such a reformulation as MILP-D often results in a binary mathematical program with exponential or potentially infinite many constraints. However, it may be tractable in the presence of special structure. If inequalities (58) preserve the totally unimodularity of the network matrix in (3)-(4), and are separable in polynomial time, then MILP-D is solvable in polynomial time in D R via the Ellipsoid method (Grötschel et al. 2012). Otherwise, computational approaches often rely on decomposition methods that iteratively add variables and constraints to the MILP.

State-augmenting Algorithm
We propose a state-augmenting approach where each iteration is a traditional combinatorial constrained shortest-path problem and is amenable, e.g., to combinatorial methods such as the pulse method (Cabrera et al. 2020). Our solution is based on recursive model DP-D as applied to RO.
We recall that, if solved via the shortest-path representation LP-D, an optimal solution corresponds to a sequence (ṡ 1 , a 1 ), . . . (ṡ n , a n ) where (a 1 , . . . , a n ) is an r − t path in D R andṡ i ∈ S ta i is a reachable state at node t a i , i ∈ {1, . . . , n}. We formalize two properties of the DP model based on this property, which build on the intuition that states conceptually play the role of constraints in recursive models. In particular, we show that only a finite set of realizations are required in the state representation to solve the Bellman equations exactly, providing a "network counterpart" of similar results in infinite-dimensional linear programs.
Proposition 9. Let z(∆ ) be an optimal solution of DP-D when the state vectorṡ is written only with respect to a subset of realizations ∆ ⊆ ∆, adjusting the dimension of d appropriately. Then, (a) z(∆ ) ≤ z(∆ ) for any ∆ ⊆ ∆ ⊆ ∆.
Proof of Proposition 9. Property (a) follows because DP-D is equivalent to RO for any ∆ ⊆ ∆ due to Proposition 8 (correctness of the reformulation) and Proposition 2 (correctness of DP-D).
Thus, since RO written in terms of ∆ has less constraints than the same model written in terms of ∆ , we must have z(∆ ) ≤ z(∆ ).
For property (b), let P be the finite set of r − t paths in D R . Partition P into subsets P F and P I = P \ P F , where p ∈ P F if and only if its associated solution x p is feasible to RO. Thus, for eacĥ p ∈ P I , there exists at least one scenario δp ∈ ∆ where some inequality of A(δp) ≤ b(δp) is violated.
We build the set which is also finite and must be such that z(∆ * ) = z(∆), since only paths in P F are feasible to CSP-D which, by definition, are also feasible to RO.
Proposition 9 and its proof suggest a separation algorithm that identifies a sequence of uncertainty sets ∆ 1 ⊂ ∆ 2 ⊂ · · · ⊂ ∆ k = ∆ * providing lower bounds to RO in each iteration, identifying scenarios to compose ∆ * . We describe it in Algorithm 1, which resembles a Benders decomposition applied to a recursive model, i.e., the state variables of LP-D encode the Benders cuts and the master problem is solved in a combinatorial way D via a constrained shortest-path algorithm.
The algorithm starts with ∆ 1 = ∅ and finds an optimal r − t path on D R . Notice that such path can be obtained by any shortest-path algorithm because the state space is initially empty. Next, the separation oracle S (x * ) is invoked to verify if there exists a realization δ for which inequalities (55) are violated. If that is the case, the state vectorṡ is augmented and the problem is resolved with a constrained shortest path algorithm. Otherwise, the algorithm terminates.
In each iteration of the algorithm, f (x * ) provides a lower bound on the optimal value of RO according to Proposition 9. Once S (x * ) = ∅, we obtain a feasible solution x * to RO and, thus, an optimal solution to the problem.
We provide a formal proof of the convergence of the algorithm below. Specifically, the worst-case complexity of the state-augmenting algorithm is a function of both (i) the number of paths in D R , (ii) the worst-case time complexity of the separation oracle S (x), and (iii) the worst-case complexity of the solution approach to LP-D to be applied. We observe that, computationally, the choice of each δ to include plays a key role in numerical performance, which we discuss in our case study of §6.3.
Proposition 10. Algorithm 1 terminates in a finite number of iterations. Moreover, the number of calls to the separation oracle S (·) is bounded above by the number of infeasible paths of D R with respect to (58).

Case Study: Robust Traveling Salesperson Problem with Time Windows
We investigate the separation algorithm on a last-mile delivery problem with uncertain service times and model it as a robust traveling salesperson problem with time windows (RTSPTW). Let where l j , u j ≥ 0 are lower and upper bounds on the service time of vertex j, respectively, and b is an uncertainty budget that controls how risk-averse the decision maker is . More precisely, lower values of b indicate that scenarios for which all vertices of V will have high service times are unlikely and do not need to be hedged against. Converserly, larger values of b imply that the decision maker is more conversative in terms of the worst case; in particular, b → +∞ represents a classical interval uncertainty set.
The objective of the RTSPTW is to find a minimum-cost route (i.e., a Hamiltonian path) starting at vertex 0 and ending at vertex n that observes time-window constraints with respect to the edge travel times. We formalize it as the MILP j:(j,i)∈E In the formulation above, the binary variable x ij denotes if the path includes edge (i, j) ∈ E and M = d n is a valid upper bound on the end time of any service. The variable w δ j is the time that vertex j ∈ V is reached when the realization of the service times is δ ∈ ∆. Inequalities (60)- (61) ensure that each vertex is visited once. Inequalities (62)-(63) represent the robust constraint and impose that time windows are observed in each realization δ ∈ ∆. Note that, since ∆ = ∅, they also rule out cycles in a feasible solution.
The main challenge in MILP-TSP is the requirement that routes must remain feasible for all realizations of the uncertainty set. In particular, solving MILP-TSP directly is not viable because of the exponential number of variables and constraints. Thus, existing routing solutions are not trivially applicable to this model, as analogously observed in related robust variants of the traveling salesperson problem (e.g., Montemanni et al. 2007, Bartolini et al. 2021).
Moreover, for simplicity of exposition, we consider that D represents directly a traditional DP formulation of the traveling salesperson problem (TSP), i.e., with each node u ∈ N we associate the state (S u , L u ) where S u contains the set of visited vertices on the paths ending at u, and L u ∈ V gives the last vertex visited in all such paths. In particular, (S u , L r ) = ({0}, 0) and L t = (V, n).
We refer to Example 4 for an instance with five vertices and to Cire and Van Hoeve (2013) for additional details and compression techniques for D.
Example 4. Figure 5 depicts D for a TSP with vertex set V = {0, . . . , 4} where every path must start at 0 and end at 4. The labels on each arc denote ν a . For example, in the left-most node of the third layer, the state (S u , L u ) = ({0, 1, 2}, 2) indicates that vertices in {0, 1, 2} were visited in all paths ending at that node, and that the last vertex in such paths is 2. Notice that the only outgoing arc a has value ν a = 3, since (i) this is the only vertex not yet visited in such path, and (ii) paths cannot finish at vertex 4.

Figure 5
A decision diagram for a TSP problem with n = 4 cities for Example 4.
The side constraints are an arc-based representation of (62). While it is possible to write an inequality system (2) in terms of arc variables y by incorporating auxiliary variables, we leverage the DP perspective to obtain a more compact state representation. Specifically, for each scenario δ ∈ ∆ and node u ∈ N , let e u (δ) be the earliest time any route finishes serving vertex L u when considering only r − t paths that include u. Formally, e v (δ) = max r Lv , min The first equality states that the earliest time we finish serving vertex 0 is zero by definition. In the second equality, the "max" term corresponds to the earliest time to arrive at vertex L v for a node v, which considers its release time and all potential preceding travel time from a predecessor vertex L u . It follows that an r − t path is feasible if and only if and the vector ( e(δ) ) δ∈∆ are the node states augmented by Algorithm 1 per realization δ.
Finally, we require a separation oracle S (x) to identify a realization violated by x. We propose to pick the "most violated" scenario in terms of the number of nodes for which the deadline is not observed (see Appendix H).
6.3.2. Numerical Study. We compare two methodologies for the RTSPTW to evaluate the benefits of the combinatorial perspective provided by DDs. The first is akin to classical approaches to RO and consists of applying Algorithm 1 using MILP-TSP when deriving a new candidate solution (Step 5); we denote such a procedure by IP. The second approach is Algorithm 1 reformulated via decision diagrams, as developed in this section. We use a simple implementation of the pulse algorithm to solve each iteration of the constraint shortest-path problems over D (see description in Appendix G). We denote our approach by DD-RO. As before, we focus on graphical description of the results, with all associated tables and detail included in Appendix I. The MILP models were solved in ILOG CPLEX 20.1. Source code is available upon request. Finally, we remark that both IP and DD-RO can be enhanced using other specialized TSPTW and constrained shortest-path formulations; here we focus on their basis comparison with respect to Algorithm 1.
Our testbed consists of modified instances for the TSPTW taken from Dumas et al. (1995).
In  Dashed lines (blue) mark the time-limit coordinate at 3,600 seconds. In total, IP and DD-RO solve average 241.9 with a standard deviation of 541.3 seconds, while the runtime for DD-RO for the same instances is 22.1 seconds with a standard deviation of 82.1 seconds. Figure 6a suggests that the complexity of instances is highly affected by time-window width for both methods. In particular, small values of w lead to smaller decision diagrams, and hence faster pulse runtimes. We observe an analogous pattern in Figure 6b, where larger values of n reflect more difficult instances. On average, runtimes for DD-RO were at least 75 times faster than IP, which underestimates the real value since many instances were not solved to optimality by IP.
The reason follows from the time per iteration between using pulse and MILP in Algorithm 1.
While both theoretically have the same number of iterations in such instance, pulse is a purely combinatorial approach over D and solves the constrained shortest-path problem, on average, in 40 seconds per iteration. For IP, the corresponding model MILP-TSP (with fewer scenarios) is challenging to solve, requiring 234 seconds per iteration on average.

Conclusion
We study reformulations of discrete bilevel optimization problems and robust optimization as We propose two approaches to solve the underlying CSP based on each setting. For the bilevel case, we leverage polyhedral structure to rewrite the original problem as a mixed-integer linear programming model, where we convexify non-linear terms connecting the master's and follower's variables by exploiting duality over the network structure. For the robust case, we proposed a stateaugmenting algorithm that iteratively add labels (i.e., resources) to existing dynamic programming approaches to CSPs, which only rely on the existence of a separation oracle that identifies violated realizations of a solution.
The methods were tested on a competitive project selection problem, where the follower and the leader must satisfy knapsack constraints, and on a robust variant of a traveling salesperson problem with time windows. In both settings, numerical results suggested noticeable improvements in solution time of the CSP-based techniques with respect to state-of-the-art methods on the tested instances, often by orders of magnitude.
Appendix A: Proof of Proposition 1 We reduce the knapsack problem to a formulation of CSP-D with the required size. Consider any for a, c ∈ R N and b ∈ R * + . Consider now a one-width decision diagram having all n permutations of {0, 1}, that is, D has one node per layer and every node except t has two ougoing arcs, one per each value assignment ν a ∈ {0, 1}. For example, the figure below depicts D for n = 4: Since D encodes the binary set {0, 1} n , it suffices to define side inequalities (2) that replicate the knapsack constraint, i.e., we fix m = 1, coefficients G 1a = c τa ν a for all a ∈ A, and d = b. Finally, for the objective, we have the lengths l a = −p τa ν a for all a ∈ A, assumed negative since CSP-D is a "min" problem (or, equivalently, we could also change the objective sense to "max"). Thus, there is a one-to-one correspondance between a knapsack instance and a solution to CSP-D, and their optimal solution values match.

Appendix B: Proof of Proposition 3
Consider the r − t arc-specified path p = (a 1 , a 2 , . . . , a n ) in D F extracted from an arbitrary feasible solution y of CSP-D. The associated solution vector x p = (ν a 1 , ν a 2 , . . . , ν an ) satisfies A F x p ≤ b F by construction of D F . Further, because of (16), we can only have x p j = ν a j = 1 for some j ∈ {1, . . . , n} if x L j = 0, that is, x p ≤ 1 − x L . Thus, since the solution values and path lengths match by construction of D F , an optimal solution of the corresponding CSP-D satisfies (13), and the converse holds analogously.
Appendix C: Proof of Proposition 7 From Proposition 5 we have that valid M -values should be large enough to provide an upper bound for the γ-variables. Consider any feasible solution to MILP-DB given by (x L ,x F ,ŷ,π,γ) and note that for all arcs a ∈ A such thatx L τa = 1 and ν a = 1, the corresponding dual γ-variables could take any value in the rangeγ a ∈ [l a +π ha −π ta , ∞], preserving feasibility. Assume without loss of generality that γ a = l a +π ha −π ta , ∀a ∈ A : x L τa = 1 and ν a = 1, that is,γ a is the shadow price corresponding to constraints (31) and captures the change in the objective function associated with adding one unit of capacity to arc a. The values forπ ha andπ ta are computed as the length of the best completion path from nodes h a and t a to t, if a completion path exists, or some arbitrary value otherwise. The relatively complete recourse property ensures that there exists at least one path from r to t and as a result Combining (68) and (69) we obtain that completing the proof.
Appendix D: Proof of Proposition 8 By construction of D R , there is a one-to-one correspondence between paths of D R and solutions of X . Consider the r − t arc-specified path p = (a 1 , a 2 , . . . , a n ) in D R extracted from an arbitrary feasible solution y of CSP-D, and let x p = (ν a 1 , ν a 2 , . . . , ν an ) ∈ X be its associated solution. Then, for any δ ∈ ∆ and i ∈ {1, . . . , r}, That is, x p is feasible to RO. Conversely, any feasible solution to RO can be shown feasible to CSP-D by applying the same steps, noticing that the objective functions are equivalent.

Appendix E: Proof of Proposition 10
At any iteration t, the algorithm either (i) stops, in which case the solution is optimal since ∆ t+1 ⊂ ∆; or (ii) x * is infeasible due to some scenario δ. Because δ belongs to the future sets ∆ t+2 , ∆ t+3 , . . . , we must have that x * is infeasible to LP-D at any future iterations, i.e., its associated path is never selected as an optimal r − t path. Thus, since the number of paths of D R is finite, the algorithm must terminate after all infeasible paths have been exhausted.

Appendix F: Decision Diagram Modeling
In this section, we briefly discuss the two general steps of modeling a problem as a discrete optimization problem DO; additional details and examples are provided in Bergman et al. (2016).
Specifically, the first step rewrites the problem as the recursive model that is, s j encodes the knapsack load at the j-th stage, T i (s j , x j ) = s j + a j x j ,f (s j , x j ) = c j x j , and Given the recursive formulation above, we write the initial D as the underlying state-transition graph, i.e., the j-th layer N j includes one node per possible value of s j , and there exists an arc with label x j between two nodes s j , s j+1 if and only if s j+1 = s j + a j x j and x j ∈ X (s j ). The length of such an arc is set tof (s i , x i ). For example, consider the knapsack instance from Example 1: s.t. 7x 1 + 5x 2 + 4x 3 + x 4 ≤ 8, The resulting (non-reduced) D is as follows, where node labels correspond to the states s j and arc labels tof (s i , x i ). We also represent all possible values of s n+1 into a single terminal node t with no loss of information. For exposition purposes, we also omit the labels on the last arcs, which are all equivalent (0 for dashed arcs, and 1 for solid arcs). Finally, the second step of the procedure is known as reduction, where we compress the graph by eliminating redundant nodes. For example, all nodes in the fourth layer can be compressed into a single node, since they have equivalent paths to t. More precisely, two nodes u, u ∈ N j are redundant if their outgoing arcs are equivalent, that is, for each arc a outgoing from u, there exists another arc a with the same head node, length, and value; and vice-versa for node u . The reduction procedure is a bottom-up traversal that merges nodes if their outgoing arcs are equivalent, i.e., 1. For each layer j = n, n − 1, n − 2, . . . , 2: (a) If there are two redundant nodes u, u ∈ N j , then merge them into a single node, redirecting the corresponding incoming arcs. Figure 1(a) depictes the resulting graph after applying the reduction mechanism above.

Appendix G: The Pulse Algorithm
The pulse algorithm proposed originally by Lozano and Medaglia (2013), is a recursive enumerative search based on the idea of propagating a pulse through a directed network, starting at the root node. Pulses represent partial paths and propagate through the outgoing arcs from each node, storing information about the partial path being explored. Once a pulse corresponding to a feasible path reaches the end node, the algorithm updates a primal bound on the objective function, stops the pulse propagation, and backtracks to continue the recursive search resulting in a pure depthfirst exploration. Opposite to labelling algorithms, pulses are not stored in memory but passed as arguments in the recursive search function. If nothing prevents pulses from propagating, the pulse algorithm enumerates all feasible paths from root node to end node ensuring that an optimal path is found.
To avoid a complete enumeration of the solution space, the pulse algorithm employs a set of pruning strategies proposed to prune pulses corresponding to partial paths, without cutting off an optimal solution. The core strategies in the pulse algorithm are pruning by infeasibility, bounds, and dominance. The infeasibility pruning strategy uses a dual bound on the minimum resource required to reach the end node from any given intermediate node and discards partial paths that cannot be completed into a feasible complete path. The bounds pruning strategy uses dual bounds on the minimum-cost (maximum-profit) paths from any intermediate node to the end node and discards partial paths that cannot be part of an optimal solution. The dominance pruning strategy uses a limited number of labels to store information on partial paths explored during the search and discard partial paths based on traditional dominance relationships. In contrast to labeling algorithms, the pulse algorithm does not rely on extending every non-dominated label for correctness; instead, this is guaranteed by properly truncating the recursive search. Hence, even if no labels were used at any of the nodes, the algorithm remains correct.
The pulse algorithm has been extended for the elementary shortest path problem with resource constraints , the biobjective shortest path problem , the weight constrained shortest path problem with replenishment (Bolívar et al. 2014), the orienteering w j = max{r j , w i + δ i + t ij }, ∀(i, j) ∈ E : x ij = 1, δ ∈ ∆, v j ∈ {0, 1}, ∀j ∈ V.
In our separation model, v j = 1 if and only if the time window of vertex j is violated, j ∈ V. The equalities (71)-(72) specify the arrival times w that follow from x. The inequalities (73) enforce the definition of v j considering a sufficiently small (e.g., the machine precision). Finally, (74) specify the domain of variables. We note that SEP-TSP can be solved by MILP solvers after linearizing the "max" in (72) with additional binary variables.
Appendix I: Additional Tables   I.1. Competitive Project Scheduling Table 1 presents extended results for our experiments with competitive scheduling instances. The first column indicates the number of items. The second column presents the budget tightness parameter. Columns 3 to 8 present for each algorithm the average computational time in seconds calculated among the 5 instances of the same configuration, the number of instances solved within the time limit, and the average optimality gap computed over the instances not solved within the time limit. We use a dash in the "Gap" column to indicate that all instances in the row are solved to optimality. For instances not solved within the time limit we record a computational time of 3600 seconds. Columns 9 and 10 present for the DD-based approach the average number of nodes and arcs in the resulting DD.  Table 2 presents additional results obtained by fixing the number of projects to n = 30 and exploring different parameters for the distribution of the projects cost as well as a wider range of the budget tightness parameter.  Table 3 presents additional results for our experiments with robust traveling salesman instances.
The first column indicates the number of nodes. The second column presents the time windows width parameter and the third column shows the budget for the uncertainty set. Columns 4 to 9 present for each algorithm the average computational time in seconds calculated among the 5 instances of the same configuration, the number of instances solved within the time limit, and the average optimality gap. Table 3 Comparing IP and DD-RO on a set of robust traveling salesman problem instances