# Generic Machine-Learning-Augmented Beam Search for Resource-Constrained Shortest Path reformulations of Combinatorial Optimization Problems

Fulin Yan, François Clautiaux, Aurélien Froger, Boris Albar

► To cite this version:

Fulin Yan, François Clautiaux, Aurélien Froger, Boris Albar. Generic Machine-Learning-Augmented Beam Search for Resource-Constrained Shortest Path reformulations of Combinatorial Optimization Problems. Computers and Operations Research, 2025, 10.1016/j.cor.2025.107339. hal-05176251v3

HAL Id: hal-05176251

https://inria.hal.science/hal-05176251v3

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Generic Machine-Learning-Augmented Beam Search for
Resource-Constrained Shortest Path Reformulations of
Combinatorial Optimization Problems

Fulin Yan¹,², François Clautiaux¹, Aurélien Froger¹, and Boris Albar²

¹Université de Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, Talence, F-33400, France

²Centre Aquitain des Technologies de l'Information et Electroniques, Talence, F-33400, France

**Abstract**

In this work, we propose a generic heuristic for the resource-constrained shortest path problem derived from dynamic programming reformulations of hard combinatorial optimization problems. The approach is a machine-learning (ML)-augmented beam search, where an ML model serves as one of the scoring functions to select candidate paths for expansion, complementing lower and upper bound estimates. Our method offers several key advantages. First, the path features used by our model are generic and only derived from the reformulation, without relying on any additional problem-specific information beyond the graph structure and resource constraints. Second, we manually designed and aggregated features to obtain vectors of fixed length, enabling us to train the model on small to medium-sized instances and apply it to much larger instances. We evaluate our algorithm on two benchmark problems: the Single Machine Total Weighted Tardiness Problem and the Temporal Knapsack Problem, each under two settings: with and without access to optimal Lagrangian multipliers. Our numerical experiments show that integration of ML into an anytime beam search enhances its solution quality in most cases, with at most minor performance degradation in the other cases. They suggest that one should systematically incorporate ML into the approach.

**Keywords:** resource-constrained shortest path reformulation, beam search, machine learning, dynamic programming

# 1 Introduction

Dynamic programming (DP) offers a powerful framework for tackling many combinatorial optimization (CO) problems by transforming them into shortest path problems within state-transition graphs. Vertices correspond to reachable states, and arcs represent transitions between states. When the resulting graph is of manageable size, typically polynomial in the size of the original problem, the shortest path can be computed efficiently, as the problem is solvable in linear time on directed acyclic graphs (DAG). In such cases, DP yields efficient solution algorithms for the original CO problem.

For NP-hard problems solvable by DP, explicitly constructing the state-transition graph is not an option due to its exponential or pseudo-polynomial size. To overcome this limitation, state-space relaxation (Christofides, Mingozzi, and Toth 1981) generates relaxed formulations with reduced state spaces. This approach yields a smaller graph that includes many infeasible paths. To ensure feasibility, additional constraints—such as capacity limits or pairwise disjunctions—must be incorporated.
Adding these constraints leads to resource-constrained shortest path problems in DAG (RCSPP-DAG), known to be NP-hard in general. Several methods have been proposed to solve them, ranging from exact methods (e.g., Ibaraki 1987; Boland, Dethridge, and Dumitrescu 2006; Lozano and Medaglia 2013; Thomas, Calogiuri, and Hewitt 2019; Cabrera et al. 2020; Himmich, El Hallaoui, and Soumis 2024) to heuristics (e.g., Avella, Boccia, and Sforza 2002; Cabrera et al. 2020; Kuroiwa and Beck 2023). We refer to Pugliese and Guerriero 2013 for the last survey on exact methods for RCSPPs.

Several exact methods (e.g., Ibaraki 1987; Dumitrescu and Boland 2003) start by building a small graph using state space relaxation, and then solve the problem obtained by Lagrangian relaxation of the resource constraints. This relaxation allows for the elimination of many arcs or partial paths that cannot belong to an optimal solution. The efficiency of these methods depends not only on the quality of the relaxation but also on the initial primal bound. In this paper, we focus on computing an initial primal bound for RCSPP-DAG.

Successful heuristics for RCSPPs employ truncated enumeration strategies. They begin at the source vertex and iteratively extend paths until they reach the sink vertex. They differ primarily from the specific criteria employed to select the paths that are retained at each iteration. Among these heuristics, beam search (BS) is particularly well-suited to problems whose solutions involve a sequence of decisions, as is the case for RCSPP-DAG. It has been successfully applied to compute high-quality solutions for diverse optimization problems (Akeb, Hifi, and M'Hallah 2009; Ghirardi and Potts 2005; Bennell and Song 2010; Blum and Miralles 2011; Libralesso et al. 2022). It is parameterized by a *beam size*, which defines the maximum number of paths maintained after each extension step. Specifically, it extends all current paths along all feasible outgoing arcs. If the number of resulting paths exceeds the beam size, a selection process discards the less promising ones. Anytime beam search (ABS) performs iterative runs of BS until the time limit is reached. At each iteration, the beam size is doubled to enlarge the search space and improve solution quality progressively. Path selection is therefore a key ingredient of BS. Most methods from the literature employ a scoring function based on either a lower bound or an upper bound for the cost of completing the paths.

For any given CO problem that admits a RCSPP-DAG reformulation, we aim to train an ML-based scoring function and integrate it within ABS to complement guidance provided by lower and upper bounds, thereby improving solution quality. In our setting, the lower bound is derived from Lagrangian relaxation of the resource constraints, while a greedy rollout algorithm computes the upper bound. The ML-based function scores and selects path extensions based solely on path information and the input data of an RCSPP-DAG instance, namely the graph structure and the resource constraints. We face two main challenges: (1) handling the varying size and structure of input graphs, and (2) designing an effective learning pipeline, including data collection, feature engineering, and model training.

RCSPP-DAG instances vary in the number of arcs, vertices, and resources, even for the same CO problem. A first challenge is designing an input representation and ML architecture that generalizes across such variable-sized instances. Unlike prior works that extract features directly from the original problem formulation (Huber and Raidl 2022; Narita, Kuroiwa, and Beck 2025), our approach must derive all features solely from the RCSPP-DAG. As a result, automatic feature encoding by a neural network (NN) is not directly applicable unless a suitable architecture is tailored to support input representations with varying structure and size. In particular, the varying number of resource constraints across instances, which leads to variable-sized arc features, combined with the complexity of encoding constraint parameters such as right-hand sides and signs, prevents the direct application of graph neural networks (GNNs). Nested architectures, such as Deep Sets (Zaheer et al. 2017) and the Set Transformer (Lee et al. 2019), can handle variable-sized sets in a permutation-invariant manner, assuming fixed-size vectors represent the elements within these sets. The heterogeneous nature of the input data in RCSPP-DAG instances complicates encoding and restricts the applicability of such architectures to our problem. Reflecting on these observations, we manually implement the
computation and aggregation of generic features into fixed-size vectors, independent of instance size.
This aggregation ensures compatibility across instances of varying dimensions and enables training
a Feed-Forward NN (FNN).

A second challenge lies in training the ML model. A key difficulty is ensuring that the solution to the learning problem is non-trivial—that is, it should not overly emphasize the dual bound score to the point where the model simply replicates it. The ML model should also capture patterns that generalize to unseen instances. Achieving both of these goals critically depends on the quality of the training data. We design a data collection process that generates partial paths along with estimates of their completion costs. These paths should be diverse, ideally include optimal solutions, and have cost estimates close to their true values. However, accurately estimating these costs is difficult, as completing a path optimally is as hard as solving the original problem. To address this, we propose techniques to increase path diversity and heuristics to improve cost estimation accuracy.

To the best of our knowledge, no prior work has integrated ML into a tree search algorithm while addressing the challenges of our input representation and its combination with other scoring functions. This specificity is a significant difference from prior works (Huber and Raidl 2022; Ettrich, Huber, and Raidl 2023; Narita, Kuroiwa, and Beck 2025) that have focused on fully using an ML-scoring function.

The document is organized as follows. Section 2 reviews related work from the literature. Section 3 describes our BS with Lagrangian filtering and multiple scoring functions, while Section 4 presents our ML pipeline designed to learn a scoring function for integration into BS. Section 5 presents the two use cases on which we tested our ML-augmented BS, and Section 6 reports the computational experiments conducted on these two test cases.

2 Literature review

In this section, we review the literature most relevant to our work. The RCSPP-DAG can be viewed as the problem of identifying an optimal sequence of decisions that transitions from an initial state to a final state. Accordingly, the first part of the review focuses on heuristic methods designed for this class of problems. The second part summarizes recent advances in tree search algorithms that incorporate ML techniques to solve CO problems with graph-structured inputs.

## 2.1 Heuristic algorithms for problems aiming at finding an optimal sequence of decisions

After a brief overview of heuristics specifically designed for RCSPPs, we turn our attention to tree search-based heuristics, which have proven effective in sequential decision-making problems. Throughout this section, we focus on general-purpose algorithms, excluding those specifically tailored to particular applications.

The most efficient algorithms for RCSPPs are designed to solve problems in relatively small graphs with very few resources (generally one or two, and fewer than ten) and only involving upper bound constraints on total resource consumption. These assumptions generally do not hold for the RCSPP-DAG problem, which contains a substantial number of resources and equality constraints. Most commonly, RCSPP algorithms construct a set of paths incrementally by iteratively appending new arcs to a list of paths until the sink node is reached. They generally involve path selection criteria, non-exact path dominance rules, or heuristic arc pruning (Costa, Contardo, and Desaulniers 2019). These methods also include exact algorithms yielding solutions even when terminated before convergence (Lozano and Medaglia 2013; Cabrera et al. 2020; Horváth and Kis 2016; Himmich, Ben Amor, et al. 2020; Himmich, El Hallaoui, and Soumis 2024). The aforementioned characteristics of the RCSPP-DAG (large number of resources and equality constraints) significantly weaken dominance rules and make a classical depth-first search strategy less effective in finding feasible
solutions. Consequently, existing heuristic algorithms for RCSPPs are not directly suitable for the
RCSPP-DAG problem.

Since the RCSPP-DAG can be cast as a sequential decision process, it is naturally suited for
solution using tree search algorithms. Complete tree search algorithms iteratively explore the entire
solution space, progressively identifying improved solutions until optimality is proven. In contrast,
heuristic (incomplete) tree search algorithms aim to quickly identify solutions by restricting explo-
ration to a portion of the solution space. The most popular methods for sequential decision process
are greedy search (GS), limited discrepancy search (LDS) (Harvey and Ginsberg 1995), and BS
(Lowerre 1976).

The algorithms we review rely on truncated enumerations, either over arcs or paths. They
construct complete solutions incrementally, starting from the source vertex and appending arcs at
each step. When path enumeration is truncated, tree search heuristics differ primarily in how they
expand partial solutions and in the strategies used to limit the number of candidate paths retained
at each step. GS retains only the best candidate path based on user-defined criteria. LDS extends
GS by allowing up to *D* deviations from greedy decisions, where *D* is the discrepancy parameter.
Discrepancies are computed by summing, along the path from the root, the number of times the
search deviates from the greedy choice. LDS explores candidate paths one at a time using a FIFO
stack, giving priority to those with greater depth and lower discrepancy. BS incrementally builds
paths by maintaining a list of paths, referred to as the *beam*. At each iteration, it considers all
possible arc extensions of the current paths in the *beam* and retains only the top *β* candidates,
according to a predefined scoring function. The value *β* is called the *width* of the *beam* or the *beam
size*.

Anytime heuristic search algorithms are tree search methods that quickly find a feasible solution and then continuously improve it until a time limit is reached or optimality is proven. Zhang (1998) proposed an anytime tree search algorithm based on iterative, truncated breadth-first searches. At each iteration, paths are pruned using a gradually relaxed heuristic criterion. Iterative Beam Search (IBS) (see Libralesso et al. (2022)) also performs successive iterations of beam search, starting with a beam size of $\beta = 1$ and doubling $\beta$ in each subsequent iteration. Anytime Column Progressive Search (ACPS) (Vadlamudi, Gaurav, et al. 2012) partitions the open list of paths $\mathcal{O}$ by depth, such that $\mathcal{O}_i \subset \mathcal{O}$ denotes paths at depth $i$. It operates with a depth index $i$; at each step, it expands the best $\beta$ paths from $\mathcal{O}_i$, inserts their successors into $\mathcal{O}_{i+1}$, and then increments $i$. The search restarts from paths at depth zero if any of the following conditions are met: the maximum depth is reached, a new best solution is found, or $\mathcal{O}_j = \emptyset$ for all $j \ge i$. The value of $\beta$ is initialized to one and incremented by one each time the search reaches the maximum depth. The Anytime Pack Progressive Search (APPS) algorithm (Vadlamudi, Aine, and Chakrabarti 2013) maintains two sets of paths: $\mathcal{O}_{\text{best}}$ and $\mathcal{O}_s$. At each iteration, APPS removes all paths from $\mathcal{O}_{\text{best}}$ and expands them, then inserts the top $\beta$ resulting paths back into $\mathcal{O}_{\text{best}}$, while the rest are added to $\mathcal{O}_s$. When $\mathcal{O}_{\text{best}}$ becomes empty, the top $\beta$ paths from $\mathcal{O}_s$ are transferred into $\mathcal{O}_{\text{best}}$. The search terminates once $\mathcal{O}_s$ is empty. In the implementation by Kuroiwa and Beck (2023), the search begins with $\beta = 1$, and $\beta$ is incremented by 1 each time paths are transferred from $\mathcal{O}_s$ to $\mathcal{O}_{\text{best}}$.

A key component in these algorithms is the scoring function used to guide the search. In most
cases, this function takes inspiration from the classical A\* algorithm (Hart, Nilsson, and Raphael
1968), where the score of a path *p* is $\psi(p) = c(p) + h(p)$, with *c(p)* representing the cost of *p* and
$h(p)$ representing an estimate of the cost of the best completion of path *p* into a feasible solution.
The estimation $h(p)$ can rely on either a dual or a primal evaluation. For instance, Kuroiwa and
Beck (2023) proposed several anytime tree-search based algorithms where function *h* provides a dual
bound. In contrast, Blum and Blesa (2007) compute $h(p)$ within BS by performing a depth-first
diving procedure from *p*, referred to as a rollout algorithm.
## 2.2 Integration of machine learning in tree-search algorithms

In recent years, many studies have leveraged ML to develop heuristic methods for CO (Bengio, Lodi, and Prouvost 2021; Mazyavkina et al. 2021). In the context of tree search algorithms, ML models are commonly employed as oracles to guide exploration by prioritizing promising regions of the solution space. Below, we review the most relevant works related to the integration of ML in heuristic tree search or branch-and-bound (B&B) algorithms. This part focuses on three questions: (1) the role of ML in the algorithm, (2) the input representation to the ML model, the choice of the model architecture and generalization, (3) the training paradigm of the ML model.

### Roles of machine learning

The use of ML to guide tree search algorithms has attracted significant attention, including algorithms for solving mixed-integer programs (MIPs) through B&B. These algorithms recursively decompose the search space into a tree, solving relaxations at each node until an integer solution is found and proved optimal. A major focus has been on using ML for node selection, that is, determining the next open node to branch (Yilmaz and Yorke-Smith 2021; Labassi, Chételat, and Lodi 2022). Yilmaz and Yorke-Smith (2021) propose a classifier that predicts whether the B&B algorithm should expand the current node's left child, right child, or both. Labassi, Chételat, and Lodi (2022) use an ML model to score open nodes, selecting for branching the one with the highest predicted score.

In heuristic tree search algorithms, ML is typically used to predict scores of partial solutions. Depending on the approach, the model may be used to produce either absolute scores (Choo et al. 2022; Huber and Raidl 2022; Narita, Kuroiwa, and Beck 2025) for individual partial solutions or relative rankings (Ettrich, Huber, and Raidl 2023) of partial solutions, typically those at the same search depth. Choo et al. (2022) integrate a policy network model into BS to guide the selection of partial solution extensions. At each iteration, the current policy is applied to eliminate low-probability candidate extensions. A greedy rollout procedure, conditioned on the policy, estimates a score for each remaining candidate. The top $\beta$ candidates are retained in the beam. They benchmarked their approach on the Traveling Salesman Problem (TSP), the Capacitated Vehicle Routing Problem (CVRP), and the Flexible Flow Scheduling Problem (FFSP).

Huber and Raidl (2022) and Ettrich, Huber, and Raidl (2023) use BS for the Longest Common Subsequence (LCS) problem, where the ML model directly predicts a score for each partial solution. The key difference between these two works lies in the loss function used for training: Huber and Raidl (2022) use mean squared error (MSE), whereas Ettrich, Huber, and Raidl (2023) explore several surrogate loss functions proposed by Negrinho, Gormley, and Gordon (2018). Narita, Kuroiwa, and Beck (2025) investigate the integration of ML into the anytime tree search algorithms IBS, ACPS, and APPS for solving domain-independent DP problems. In each case, a learned scoring function evaluates partial solutions and guides the search by prioritizing those to expand. They tested these algorithms on TSP, TSP with time windows (TSPTW), 0-1 Knapsack Problem, and 4-moment Portfolio Optimization Problem.

### Input data representation, model architecture, and generalization

Classical ML methods generally require fixed-size feature vectors, whereas many CO problems involve variable-sized inputs that scale with instance size. In the case of RCSPP-DAG, the number of arcs, vertices, and resource constraints grows with the instance size. Learning in this context requires input representations that capture the structural properties of the problem and model architectures that can process variable-sized data. Both elements are essential for enabling generalization to instances of sizes not encountered during training. In the following, we review how existing works address these challenges.
In B&B algorithms for MIP, ML-based node selection requires input features that effectively represent each node, which corresponds to an MIP subproblem. Yilmaz and Yorke-Smith (2021) categorize these features into three types: variable-level features (e.g., objective coefficient, fractionality), node-level features (e.g., lower bound, subtree solution estimates), and global features that describe the overall search state (e.g., best primal bound, current dual bound, integrality gap). Their FNN requires separate training for each instance size because variable-level features depend on the number of decision variables, thereby limiting generalization. To overcome this limitation, Labassi, Chételat, and Lodi (2022) represent each MIP subproblem as a bipartite graph, enabling the use of GNNs (see Wu et al. (2022)). One set of nodes corresponds to constraints, the other to variables, with edges connecting variables to constraints whenever the variable has a non-zero coefficient. Constraint nodes carry features like right-hand side values, while variable nodes include objective coefficients, bounds, and types. Edge features capture the non-zero coefficients linking variables and constraints. GNNs naturally handle varying graph sizes, thus allowing generalization across different MIP instances. However, they require that all feature vectors have fixed, instance-independent dimensions.

In heuristic tree search algorithms, we review two recent works that have studied variable instance sizes (Narita, Kuroiwa, and Beck 2025; Choo et al. 2022). They differ from previous works, which were typically problem-specific and dependent on instance size (e.g. Huber and Raidl (2022)). In both works, different input representations are employed for different problems, allowing the use of a problem-specific ML model architecture to handle variable instance sizes. For TSP and CVRP, Choo et al. (2022) represent the input as a set of fixed-size feature vectors, each corresponding to a vertex in the transportation graph. They use the transformer-based architecture proposed by Kool, Hoof, and Welling (2019), which is specifically designed to handle variable-sized input sets for routing problems. As a result, their ML model demonstrates the ability to generalize to larger instances. For FFSP, the authors represent the input as a set of matrices, each corresponding to a stage, where each matrix encodes the processing time of a job on a machine at that stage. Each matrix can be interpreted as a bipartite graph: one set of nodes corresponds to jobs, the other to machines, and edges represent the processing times. This bipartite graph structure is well-suited for processing with GNNs. Consequently, the authors adopt the model architecture proposed by Kwon et al. (2021), where a distinct GNN is applied to each stage's bipartite graph. The resulting ML model is independent of the number of jobs and machines but depends on the number of stages. The authors evaluate the generalizability of the model under a fixed number of stages and machines, and show that the model generalizes well as the number of jobs increases. Narita, Kuroiwa, and Beck (2025) also explore instance-size-independent model architectures across four combinatorial problems. For TSP and TSPTW, it exploits the inherent graph structure, where vertices represent customers and edges correspond to possible travel paths. Fixed-size node and edge features enable the use of GNNs. In the 0-1 Knapsack problem, items are encoded as variable-sized sets of fixed-dimension feature vectors (weight, profit, selection flag), which are processed using Deep Sets. Similarly, in the 4-moment portfolio optimization problem, investments are represented by fixed-size vectors (cost, expected return, variance, selection status), with the variable number of investments handled using the Set Transformer architecture. While Narita, Kuroiwa, and Beck (2025) employ instance-size-independent models, their ability to generalize to larger instances remains, to our knowledge, untested.

## Training Paradigms

The ML models used to guide tree search methods differ in the learning paradigm used (online or offline), in the architecture used, and in the features selected. Online methods are generally based on *reinforcement learning* (RL), which involves an agent interacting with an environment to iteratively improve a policy through exploration and exploitation of informative feedback. Huber and Raidl (2022) and Ettrich, Huber, and Raidl (2023) train value-based models. At each epoch,
they generate new instances and use BS, guided by the ML model, to sample promising paths. For each collected path, a nested BS computes a primal bound on the cost of completing the path. The ML model is then updated using the newly generated pairs of (path, cost). Choo et al. (2022) use policy-based RL, while Narita, Kuroiwa, and Beck (2025) use value-based and policy-based RL.

The offline learning paradigm is often employed to incorporate ML models into B&B algorithms (Labassi, Chételat, and Lodi 2022; Yilmaz and Yorke-Smith 2021). It involves two main steps: (1) node sampling and (2) supervised training on the sampled data. In the sampling phase, an MIP solver solves the training instances to optimality and also identifies the top-k solutions. A B&B algorithm, guided by a specialized node selection heuristic, is then run to collect representative nodes. Yilmaz and Yorke-Smith (2021) and Labassi, Chételat, and Lodi (2022) collect nodes along paths to an optimal solution. Yilmaz and Yorke-Smith (2021) extends this by sampling nodes along the paths to the top-k paths. To increase dataset diversity, Labassi, Chételat, and Lodi (2022) also sample nodes that slightly deviate from optimal paths.

# 3 Beam search with Lagrangian filtering and multiple scoring functions

This section presents a tailored BS for a RCSPP-DAG resulting from the reformulation of a CO problem using DP. The algorithm uses multiple composite scoring functions and graph filtering techniques based on lower and upper bounds. We begin with a formal problem definition, followed by a description of our specialized BS and ABS, a discussion of a Lagrangian filtering strategy, and an explanation of the use of various scoring functions.

## 3.1 The DPRC(P) problem: definition

Let $P$ be a combinatorial optimization problem that admits an RCSPP-DAG formulation via dynamic programming. We refer to the resulting problem as DPRC($P$) (short for dynamic program with resource constraints of $P$). It consists of a graph and a set of resource constraints. Let $\mathcal{G} = (\mathcal{V}, \mathcal{A})$ be a directed acyclic multigraph, where $\mathcal{V}$ is the set of vertices representing the states and $\mathcal{A}$ the set of arcs representing the transitions between them. Each arc $a \in \mathcal{A}$ has a tail $v^{-}(a) \in \mathcal{V}$, a head $v^{+}(a) \in \mathcal{V}$, and a cost $c_a \in \mathbb{R}$. We denote by $\mathcal{A}^{-}(v)$ and $\mathcal{A}^{+}(v)$ the sets of arcs that enter and leave the vertex $v$, respectively. The set $\mathcal{V}$ contains a designated source $s$ with no incoming arcs and a designated sink $t$ with no outgoing arcs. A path $p$ is a sequence of arcs $(a_1, \dots, a_n)$ such that $v^{-}(a_i) = v^{+}(a_{i-1})$ for $i = 2, \dots, n$. Its cost $c(p)$ equals $\sum_{a \in p} c_a$. For any path $p$, we denote by $v^{-}(p)$ and $v^{+}(p)$ the starting and ending vertices of $p$. We call $u-v$ path a path $p$ such that $v^{-}(p) = u$ and $v^{+}(p) = v$. The number of arcs in a path $p$ is referred to as its size and denoted by $|p|$. For every $u, v \in \mathcal{V}$, let $\mathcal{P}_{u,v}$ denote the set of $u-v$ paths.

The problem definition also involves resource constraints. Each resource constraint $k \in \mathcal{R}$ has a capacity denoted by $R^k$, and a consumption $r_a^k \in \mathbb{Z}$ defined on each arc $a \in \mathcal{A}$. The resource consumption of a path $r(p) \in \mathbb{Z}^{|\mathcal{R}|}$ is the cumulative resource consumption of its arcs, i.e. for a given path $p$, and a given resource $k$, $r^k(p) = \sum_{a \in p} r_a^k$. The set $\mathcal{R}$ is partitioned into two subsets: $\mathcal{R}^=$, containing equality constraints, and $\mathcal{R}^\le$ containing inequality constraints.

**Problem 1 (DPRC(P)).** Given a pair $(\mathcal{G}, \mathcal{R})$, where $\mathcal{G}$ is a DAG, and $\mathcal{R}$ a set of resource constraints, the DPRC($P$) problem consists of finding among the s-t-paths $p = (a_1, \dots, a_n)$ in $\mathcal{G}$ satisfying $r^k(p) \le R^k, \forall k \in \mathcal{R}^\le$ and $r^k(p) = R^k, \forall k \in \mathcal{R}^=$ a path that minimizes $c(p)$.

Note that the problem is parametrized by a problem $P$. We do not make any assumptions about $P$, except that it admits a DP with side constraints reformulation. Constraints are organized into families indexed by a set $\mathcal{F}$. For each family $f \in \mathcal{F}$ of resource constraints, the number of constraints in $f$ is instance-dependent. However, all constraints within $f$ share the same right-hand
side. For example, for set-partitioning constraints, each element to cover leads to a constraint of the form “= 1". Furthermore, the minimum and maximum resource consumption of any s-u path (with $u \in \mathcal{V}$) for every resource constraint in $f$ are instance-independent. We provide an MIP formulation of DPRC(P) in A.1.

## 3.2 Beam search for DPRC(P)

We now describe how we compute a solution to a DPRC(P) instance $(\mathcal{G}, \mathcal{R})$ using Beam Search. Beam Search (BS) is an iterative algorithm that incrementally extends s-v paths until reaching the sink node $t$. A central feature of BS is the use of a bounded list of paths, called the *beam*, limited to size $\beta$. At every iteration, referred to as *level*, the algorithm processes each path currently in the beam. Each path is expanded by appending all outgoing arcs from its endpoint, generating new candidate paths. Among the newly generated paths, only the $\beta$ most promising ones, according to a score, are retained for the next iteration. Throughout the search, the algorithm maintains the best feasible s-t-path found. The search ends when the beam becomes empty.

To describe the algorithm formally, we use the operator $\oplus$, which concatenates two paths to form a new path. If $p^1$ and $p^2$ are two paths of sizes $n_1$ and $n_2$ such that $v^+(p^1) = v^-(p^2)$, then $p^1 \oplus p^2 = (a_1^1, a_2^1, \dots, a_{n_1}^1, a_1^2, a_2^2, \dots, a_{n_2}^2)$. For convenience, we also introduce $p_s$ as the empty path originating at vertex $s$, i.e., $v^-(p_s) = v^+(p_s) = s$.

## 3.3 Lagrangian relaxation and graph filtering

Our customized BS performs a Lagrangian relaxation of the resource constraints from $\mathcal{R}$. Let $\lambda$ be a vector of size $|\mathcal{R}|$, in which each element $\lambda^k \in \mathbb{R}, \forall k \in \mathcal{R}^-, \lambda^k \in \mathbb{R}^+, \forall k \in \mathcal{R}^\le$ represents the Lagrangian multiplier associated with the $k$-th resource constraint.

From a given vector $\lambda$ of Lagrangian multipliers, we define the following quantities:

* the Lagrangian cost of any arc $a \in \mathcal{A}: \tilde{c}_a^\lambda = c_a + \sum_{k=1}^{|\mathcal{R}|} \lambda^k r_a^k$

* the Lagrangian cost of any path $p$ in $\mathcal{G}: \tilde{c}^\lambda(p) = -\sum_{k \in \mathcal{R}} \lambda^k R^k + \sum_{a \in p} \tilde{c}_a^\lambda$

* the minimum Lagrangian cost of any $v-t$ path with $v \in \mathcal{V}: \underline{h}^\lambda(v) = \min_{p \in \mathcal{P}_{v,t}} \{\tilde{c}^\lambda(p)\}$.

In a preliminary step, we consider two settings: (1) when the vector of the best Lagrangian multipliers, $\lambda^*$, is calculated, i.e. we set $\lambda$ to $\lambda^* \in \operatorname{argmax}_\lambda \min_{p \in \mathcal{P}_{s,t}} \tilde{c}^\lambda(p)$ and (2) when $\lambda^*$ is not calculated, i.e. we set $\lambda$ to the null vector $\mathbf{0}$. The best Lagrangian multipliers are calculated using different algorithms based on the CO problem. We will develop later the algorithms used for two CO problems in Section 6.

Given an initial upper bound $\text{UB}_{\text{initial}}$ on the value of the problem, the multipliers $\lambda$ are used to apply a reduced-cost variable-fixing method (Pessoa et al. 2010; Kramer et al. 2019). If an arc $a$ is such that $\min_{p:a \in p} \{\tilde{c}^\lambda(p)\} > \text{UB}_{\text{initial}}$, it is removed (filtered) from $\mathcal{G}$, since any path including $a$ has cost at least $\text{UB}_{\text{initial}}$. The algorithm also removes the vertices that no longer have arcs connected to them.

## 3.4 Beam search with multiple scoring functions

BS can be adapted to support multiple scoring functions applied sequentially. We define a *path selector* as a tuple $(\psi, b, e)$, where $\psi$ denotes a scoring function, $b \in \mathbb{N}$ specifies the number of paths to select, and $e \in \{0, 1\}$ indicates whether the selection uses a roulette wheel mechanism. Given a set of paths $\mathcal{P}$, the path selector $(\psi, b, e)$ proceeds as follows. If $e=0$ (deterministic path selector), the $\min\{b, |\mathcal{P}|\}$ paths with the smallest $\psi$-scores are selected. If $e=1$ (rosette wheel), $\min\{b, |\mathcal{P}|\}$
paths are selected randomly without replacement, where the probability of selecting path p from P is $\frac{\max_{p' \in P}\{\psi(p')\} - \psi(p)}{\sum_{p' \in P} \max_{p' \in P}\{\psi(p')\} - \psi(p')}.$

We now introduce Beam Search with Multiple Scoring Functions (BSM), specialized for DPRC(P). The algorithm is parameterized by a list of path selectors $\Psi = ((\psi_1, b_1, e_1), (\psi_2, b_2, e_2), \dots)$. The order of the path selectors implies that a scoring function $\psi_i$ selects paths that are complementary to those chosen by the preceding ones $\psi_{i'}$, $\forall i' \in \{1, \dots, i-1\}$.

BSM filters paths by comparing their associated Lagrangian dual bounds with the current upper bound UB. For any s-v path p, and any multiplier vector $\lambda$, we denote by $\nu^\lambda(p)$ a lower bound on the cost of any s-t path $p'$ such that $p' = p \oplus p''$ for the v-t paths $p''$ satisfying $\underline{\lambda}(v) = \tilde{c}^\lambda(p'')$. We set:

$$\nu^{\lambda}(p) = \max(\tilde{c}^{0}(p) + \underline{h}^{0}(v^{+}(p)), \tilde{c}^{\lambda}(p) + \underline{h}^{\lambda}(v^{+}(p)))$$

A path p such that $\nu^\lambda(p) \ge \text{UB}$ cannot participate in forming an improving solution.

BSM can use a greedy rollout algorithm that returns a feasible s-t path $p_{\text{rol}}$ that completes an s-v path, where $v \in V \setminus \{t\}$. By convention, $c(p_{\text{rol}}) = \infty$ if the rollout algorithm does not find a feasible completion.

Algorithm 1 outlines the general structure of BSM. The beam $\mathcal{B}$ is initialized with $p_s$ (line 1). At each step, the algorithm extends each path $p \in \mathcal{B}$ by appending every arc originating from $v^+(p)$ (line 6). Every newly created path $p'$ is stored in $\mathcal{B}'$ if $p'$ is feasible, $\nu^\lambda(p') < \text{UB}$ and $v^+(p') \neq t$ (lines 7- 12). A path $p'$ is feasible if it is part of a feasible s-t path. Determining the feasibility of $p'$ is hard in general, and thus, we use necessary conditions (see B.1). The algorithm then removes dominated paths from $\mathcal{B}'$ (line 13). A path $p_1$ dominates a path $p_2$ if any feasible s-t path containing $p_2$ would also be feasible by replacing $p_2$ with $p_1$ and would have a lower cost. We use sufficient conditions to check the dominance between paths (see B.2). If available, the rollout algorithm completes each non-dominated path $p$ into a solution to the problem (Line 17). If the cost of the path $p \oplus p_{\text{rol}}$ is better than UB, then UB is updated to $c(p \oplus p_{\text{rol}}) = \nu^\lambda(p)$, we already know an optimal completion of $p$, and we discard this path. Finally, the algorithm limits the number of paths to expand in the next depth. It first selects up to $b_1$ paths from $\mathcal{B}'$ based on the $\psi_1$-scoring function. It then selects up to $b_2$ paths that remain in $\mathcal{B}'$ based on $\psi_2$. This process continues until the final scoring function is applied or $\mathcal{B}'$ is empty (lines 22 - 27).

Our anytime beam search with multiple scoring functions (ABSM) executes BSM at each iteration, starting from a small beam size $\beta_a$. At each iteration, $\beta_a$ is doubled, following Kuroiwa and Beck (2023) and Libralesso et al. (2022), and BSM leverages the best primal bound obtained so far to prune paths using dual bounds. Early stopping is triggered once the current best solution is proven optimal. Optimality is proven if one of the following conditions holds: (1) the gap between the primal and dual bounds is zero, or (2) the current iteration corresponds to a complete path enumeration, i.e., no paths are discarded by the beam restriction.

## 4 Machine learning guided beam search for DPRC(P)

Given a CO problem P, we now explain how we leverage ML techniques to produce a scoring function $\tilde{\psi}_\omega : p \mapsto \tilde{\psi}_\omega(p) \in \mathbb{R}$ (see Algorithm 1). The scoring function $\tilde{\psi}_\omega$ is parameterized by $\omega$, a vector of learnable parameters. The goal of this scoring function is to identify paths that are complementary to those selected by two standard scoring functions. To clarify, the ML-augmented version of BSM employs a list of three path selectors, each relying on a distinct scoring function. The first scoring function (hereafter referred to as the Lagrangian dual scoring function) scores a path $p$ based on a lower bound $\nu^\lambda(p)$ on the cost of an s-t path containing $p$. The second scoring function scores a path $p$ based on the cost of an s-t path obtained from $p$ by a greedy rollout algorithm.
<caption>**Algorithm 1:** BSM((G, R), UB, p₀, Ψ) a beam search with multiple scoring functions</caption>

**Input** : Instance $(\mathcal{G}, \mathcal{R})$, initial upper bound UB, initial path $p_0$ to complete into a feasible s-t path, list of path selectors $\Psi$
**Output**: Best-known feasible s-t path ($p^*$), the paths at each level $(\{\mathcal{B}_j\}_j)$, optimal completion paths ($\mathcal{F}$)
1 &nbsp; $\mathcal{F} \leftarrow \emptyset$, $\mathcal{B}_0 \leftarrow \{p_0\}$, $p^* \leftarrow ()$, $j \leftarrow 0$;
2 &nbsp; **while** $\mathcal{B}_j \neq \emptyset$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{B}' \leftarrow \emptyset$;
4 &nbsp;&nbsp;&nbsp;&nbsp; **for** $p \in \mathcal{B}_j$ **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $a \in \mathcal{A}^+(v^+(p))$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $p' \leftarrow p \oplus (a)$;
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** IsFeasible($\mathcal{G}, \mathcal{R}, p'$) and $\nu^\lambda(p') < \text{UB}$ **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $v^+(p') = t$ **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $c(p') < \text{UB}$ **then**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor p^* \leftarrow p', \text{UB} \leftarrow c(p') \rfloor$;
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **else**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor \mathcal{B}' \leftarrow \mathcal{B}' \cup \{p'\} \rfloor$;
13 &nbsp;&nbsp;&nbsp;&nbsp; RemoveDominatedPaths($\mathcal{B}'$);
14 &nbsp;&nbsp;&nbsp;&nbsp; $j \leftarrow j + 1$;
15 &nbsp;&nbsp;&nbsp;&nbsp; **if** Rollout is available **then**
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $p \in \mathcal{B}'$ **do**
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $p_{\text{rol}} \leftarrow \text{Rollout}(p)$;
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $c(p) + c(p_{\text{rol}}) < \text{UB}$ **then**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor p^* \leftarrow p \oplus p_{\text{rol}}, \text{UB} \leftarrow c(p) + c(p_{\text{rol}}) \rfloor$;
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $c(p) + c(p_{\text{rol}}) = \nu^\lambda(p)$ **then**
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor \mathcal{B}' \leftarrow \mathcal{B}' \setminus \{p\}, \mathcal{F} \leftarrow \mathcal{F} \cup \{p_{\text{rol}}\} \rfloor$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ▷ $p_{\text{rol}}$ is an optimal completion of $p$.
22 &nbsp;&nbsp;&nbsp;&nbsp; **for** $(\psi, b, e) \in \Psi$ **do**
23 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $e = 1$ **then**
24 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{B}_\psi \leftarrow \text{select from } \mathcal{B}'$ up to $b$ paths by roulette wheel based on the scores given by function $\psi$;
25 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **else**
26 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor \mathcal{B}_\psi \leftarrow \text{select from } \mathcal{B}'$ up to $b$ paths with the smallest scores given by function $\psi$;
27 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{B}' \leftarrow \mathcal{B}' \setminus \mathcal{B}_\psi$, $\mathcal{B}_j \leftarrow \mathcal{B}_j \cup \mathcal{B}_\psi$;
28 &nbsp;&nbsp;&nbsp;&nbsp; **return** $p^*$, $\{\mathcal{B}_j\}_j$, $\mathcal{F}$;
The third scoring function is the learned function $\tilde{\psi}_{\omega}$. It is worth noting that the construction of this function depends on whether the vector of optimal Lagrangian multipliers is available, as the features differ accordingly. Hereinafter, we refer to ML-BSM and ML-ABSM as this ML-augmented version of BSM and ABSM, respectively.

We use supervised learning (SL) to train a suitable scoring function $\tilde{\psi}_{\omega}$. A significant challenge lies in constructing the training set, that is, determining a set of paths $p$ with a precise evaluation of the cost of an $s-t$-path containing $p$ (hereinafter referred to as total completion cost), if such a path exists.

Our method addresses the following two issues. Firstly, given a path $p$, determining whether it is part of a feasible $s-t$ path and, if so, determining its best total completion cost is an NP-complete problem. Secondly, it is undesirable to generate an excessive number of paths for which the Lagrangian information is already perfect (i.e., paths for which the total completion cost is equal to the lower bound). Among these paths, those for which the total completion cost, computed by a greedy rollout algorithm, equals the lower bound are the least desirable.

We construct the training pipeline shown in Figure 1. The pipeline consists of three successive steps: (1) data collection, (2) learning, and (3) testing. The data collection phase (see Section 4.1) processes instances of DPRC(P) individually to generate a set of pairs, each comprising a path $p$ and its corresponding total completion cost $y(p)$. These pairs constitute the training set for the subsequent learning phase. The learning phase (see Section 4.2) uses offline SL to predict the total completion cost of paths, yielding $\tilde{\psi}_{\omega}$. The training procedure minimizes the mean square error (MSE) between $\tilde{\psi}_{\omega}(p)$ and $y(p)$ on the collected paths. The testing phase evaluates ML-ABSM on unseen instances of DPRC(P).

<figure><img src="image_2.png" /><figcaption>Figure 1: Pipeline</figcaption></figure>

During both data collection and testing, we first apply a preprocessing step to reduce an input instance $(\mathcal{G}_0, \mathcal{R}_0)$ to an instance $(\mathcal{G}, \mathcal{R})$ containing fewer vertices and arcs. The process begins by executing BSM$((\mathcal{G}_0, \mathcal{R}_0), +\infty, p_s, ((\nu^{\lambda}, \beta_0, 0)))$ with the rollout algorithm deactivated. The initial primal bound $\text{UB}_{\text{initial}}$ is then set to the cost of the returned solution $p^*$. To ensure a fast computation of $\text{UB}_{\text{initial}}$, the beam size $\beta_0$ is chosen to be relatively small. Subsequently, the graph filtering algorithm (see Section 3.3) removes arcs and vertices that cannot be part of any solution with a cost lower than $\text{UB}_{\text{initial}}$.

## 4.1 Path sampling and score collection

The data collection phase takes as input a set $\mathcal{I}$ of instances of DPRC(P), and produces a set of (path, score) pairs. For any given instance $(\mathcal{G} = (\mathcal{V}, \mathcal{A}), \mathcal{R}) \in \mathcal{I}$, it builds a set of $s-v$ paths, where $v \in \mathcal{V} \setminus \{s\}$, ensuring that each $s-v$ path is part of at least one feasible $s-t$ path. It also computes a total completion cost for each of these paths, thereby establishing the target score of that path. Even with medium-sized instances, finding the best total completion of a path is computationally impractical. Therefore, we resort to computing upper bounds. In this phase, we employ heuristics that are more computationally intensive than those used in a greedy rollout algorithm. The objective is to learn a scoring function that produces more accurate total completion costs than this simple method.
The collection procedure consists of two main steps: path sampling and collection of total completion costs. Given a set of instances as input, the path sampling step generates a set of paths from multiple BSM runs with partial randomization. The selection of paths is based on two criteria: the potential for yielding high-quality solutions and the diversity it brings. Incorporating random paths ensures that the second criterion is met. We exclude any path p such that $c(p)+c(\text{Rollout}(p)) = \nu^{\lambda}(p)$ from the training data, since BSMF never scores it. As a byproduct of the process, the total completion costs of the generated paths originate from the rollout algorithm, backtracking, or a stitching process (as shown in Algorithm 5). The second step selects paths for the learning phase and performs BS from each path's end vertex to improve its total completion cost.

Algorithm 2 outlines the general scheme of the path sampling procedure. The algorithm takes a problem instance as input and produces a set of (path, score) pairs. For this purpose, it performs multiple runs of BSM (Algorithm 1) using three path selectors designed to introduce a diversity of paths while retaining those that could potentially lead to feasible s-t paths, thereby balancing exploration and exploitation.

To provide a precise description, we introduce the following notation. Let $h_p^{[m]}$ be the forward completion bound computed for a path p at iteration m. When we record a forward completion cost $h_p^{[m]}$, we associate this cost with all paths $p'$ such that $v^+(p) = v^+(p')$ and $r(p) = r(p')$, since they have the same completion. For any path p, and at each iteration m, let $N_p^{[m]}$ denote the number of times up to iteration m that a forward completion has been computed for all paths $p'$ sharing the same completion as p. We set $h_p^{[m]} = \infty$ if we computed no forward completion for p at iteration m.

Algorithm 2 uses the following scoring functions, each associated with a single path selector. The first scoring function, denoted $\psi_{p^*}$, ensures that we always keep in the beam the s-v paths contained in the feasible s-t path $p^*$. Specifically, the paths along $p^*$ are scored to $-\infty$, all others are scored to $\infty$. In practice, $p^*$ is the best feasible s-t path computed so far. The second scoring function, denoted $\psi_{\text{avg}}$, favors the paths that have the best total completion bound, with a larger weight for the case where the path often led to good feasible s-t paths in previous iterations. The function is precisely defined as follows:

$$ \psi_{\text{avg}}(p) = \frac{1}{1 + N_p^{[m-1]}} \left( \nu^{\lambda}(p) + \min_{i \in \{1, \dots, m-1\}} \{c(p) + h_p^{[i]}\} \times N_p^{[m-1]} \right). \quad (1) $$

The value of $\psi_{\text{avg}}(p)$ is a convex combination of the Lagrangian lower bound $\nu^{\lambda}(p)$ and the best total completion cost observed up to the previous iteration $m-1$. Since the total completion cost is always greater than or equal to $\nu^{\lambda}(p)$, we have $\psi_{\text{avg}}(p) \ge \nu^{\lambda}(p)$. The weights assigned to the two terms reflect the confidence in each estimate. If no path sharing the same completion as p has been explored yet (i.e., $N_p^{[m-1]} = 0$), $\psi_{\text{avg}}(p)$ reduces to the Lagrangian lower bound $\nu^{\lambda}(p)$, which is the minimal value $\psi_{\text{avg}}(p)$ can take, encouraging the beam search to explore it. As $N_p^{[m-1]}$ increases, the weight of the total completion cost grows, indicating that greater exploration of paths with the same completions leads to higher confidence in the total completion cost. Early large total completion costs can prevent certain paths from being explored. To alleviate this issue, we introduce a third scoring function. This scoring function, denoted $\psi_{\text{rand}}$, aims at diversity by selecting purely random paths, following a uniform distribution $U(0, 1)$.

Algorithm 2 begins with an initial run of BSM (line 4), starting from the empty path $p_s$ at the source node s. In this step, we use the Lagrangian dual bound $\nu^{\lambda}(p)$ as the sole scoring function, and paths are selected deterministically. Throughout the iterations, the algorithm populates a set of paths (set $\mathcal{D}$ in Algorithm 2) with their current minimum total completion cost (line 10), excluding any path for which the completion given by the greedy rollout algorithm is proven optimal by the Lagrangian dual bound. We repeatedly execute BSM using three different path selectors (line 6): the first selects the best path according to $\psi_{p^*}$, the second uses a roulette-wheel mechanism based on scores from $\psi_{\text{avg}}$, and the third selects paths based on $\psi_{\text{rand}}$. After running BS, the algorithm computes the forward completion cost of each collected path (line 7). We compute this cost via a
call to Algorithm 5 (detailed in C), which begins with the last beam $\mathcal{B}_J$ and proceeds backwards, down to $\mathcal{B}_1$. If more than one path can be appended to form a feasible s-t path, the algorithm selects the one with the minimum cost. At the end of each iteration, Algorithm 2 updates $\psi_{\text{avg}}$ using Equation (1), and adjusts $\beta_{\text{avg}}$ based on the number of iterations $m$ (see line 11). The value of $\beta_{\text{avg}}$ is gradually decreased as $m$ increases to introduce greater diversity among the sampled paths. The objective is twofold: to improve the forward completion costs of the already sampled paths and to explore different paths in subsequent iterations. The sampling procedure stops when it reaches the maximum number of iterations $M$. It returns a set of paths along with their minimum total completion costs. Note that it does not include s-t paths, which are never evaluated during BS.

<caption>**Algorithm 2: Path sampling**</caption>

**Input**: $(\mathcal{G}, \mathcal{R})$, $\text{UB}_{\text{initial}}, \beta, \beta_{\text{avg}}$
**Output**: a set of paths with their scores
1 &nbsp; $\text{UB} \leftarrow \text{UB}_{\text{initial}}, \mathcal{D} \leftarrow \emptyset$;
2 &nbsp; for $m \in \{0, \dots, M\}$ do
3 &nbsp;&nbsp;&nbsp;&nbsp; if $m=0$ then
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor p^*, \{\mathcal{B}_j\}_j, \mathcal{F} \rfloor \leftarrow \text{BSM}((\mathcal{G}, \mathcal{R}), \text{UB}, p_s, ((\nu^\lambda, \beta, 0)))$;
5 &nbsp;&nbsp;&nbsp;&nbsp; else
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor p^*, \{\mathcal{B}_j\}_j, \mathcal{F} \rfloor \leftarrow \text{BSM}((\mathcal{G}, \mathcal{R}), \text{UB}, p_s, ((\psi_{p^*}, 1, 0), (\psi_{\text{avg}}, \beta_{\text{avg}}, 1), (\psi_{\text{rand}}, \beta - \beta_{\text{avg}} - 1, 0)))$;
7 &nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{H} \leftarrow \text{Backtrack}((\mathcal{G}, \mathcal{R}), \{\mathcal{B}_j\}_j)$;
8 &nbsp;&nbsp;&nbsp;&nbsp; for $(p, h_p^{[m]}) \in \mathcal{H}$ do
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if $p \notin \mathcal{F}$ then
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor \text{Update } \mathcal{D} \text{ with the pair } (p, c(p) + h_p^{[m]})$;
11 &nbsp;&nbsp;&nbsp;&nbsp; Update $\psi_{\text{avg}}$ by Equation (1) and ajust $\beta_{\text{avg}}$ by $m$;
12 &nbsp; return $\mathcal{D}$

$\triangleright$ Algorithm 5

We run the path sampling algorithm on several problem instances to create our training set. To control the size of the training set, for each instance $(\mathcal{G}, \mathcal{R})$, we uniformly randomly select at most $N_{\text{enhance}}$ paths generated by the path sampling algorithm (i.e., from $\mathcal{D}$ in Algorithm 2) on each instance. For each pair $(p, q_p)$ of a path $p$ and its total completion cost $q_p$, we run $\text{BSM}((\mathcal{G}, \mathcal{R}), q_p, p, ((\nu^\lambda, \beta_{\text{enhance}}, 0))$ to obtain a better completion cost. Let $p_{\text{compl}}$ be the computed path. We set the score of path $p$, denoted by $y(p)$, to the minimum between $c(p) + c(p_{\text{compl}})$ and $q_p$. The score of a path corresponds to the minimum cost of a solution computed either by the greedy rollout algorithm, the backtracking and stitching process, or the above-mentioned BS from $v^+(p)$. This combined approach for computing a completion to a path is henceforth referred to as the *enhanced rollout algorithm*.

## 4.2 Learning a total completion cost of a path

We now describe our SL approach for predicting the total completion cost of a path. For clarity and precision, we denote $\ell$ as a instance-path pair $\ell = ((\mathcal{G}(\ell), \mathcal{R}(\ell)), p(\ell))$, where $(\mathcal{G}(\ell), \mathcal{R}(\ell))$ denotes an instance of DPRC(P) and $p(\ell)$ is a path in graph $\mathcal{G}(\ell)$. The training dataset consists of a set of instance-path pairs, each associated with a score. The score $y(p(\ell))$ assigned to a pair $\ell$ represents the *total completion cost* for $p(\ell)$ collected during the previous phase.

Our objective is to learn a scoring function $\tilde{\psi}_{\omega}$ that can be applied across different instances of DPRC(P). This function is parameterized by a vector $\omega$ whose values are determined through training. As the size of each instance - defined by $|\mathcal{V}|$, $|\mathcal{A}|$, and $|\mathcal{R}|$ - varies, the amount of information encoded in a instance-path pair depends on the instance. To use a Feed-Forward Neural Network (FNN), we create a fixed-size representation $\phi(\ell) \in \mathbb{R}^d$ of each instance-path pair $\ell$. We represent
the FNN as a function $\tilde{\psi}(\cdot, \omega)$ parameterized by $\omega$, which maps any feature vector $\mathbf{f} \in \mathbb{R}^d$ to a real-valued score: $\tilde{\psi}(\mathbf{f}, \omega) \in \mathbb{R}$. We formally define the scoring function for a path associated with an instance-path pair $\ell$ as follows:

$$ \tilde{\psi}_{\omega}(p(\ell)) = \tilde{\psi}(\phi(\ell), \omega) \qquad (2) $$

The learning phase consists of two main steps: feature extraction and training. First, the feature extraction step constructs a mapping $\phi$ from an instance-path pair $\ell$ to a vector $\phi(\ell)$, in which each element represents a feature of $\ell$. Second, the training step uses SL to train a neural network to predict the score of an instance-path pair based on its feature vector.

### 4.2.1 Feature extraction

The feature extraction step involves extracting a fixed-size representation of each instance-path pair $\ell$. We first describe the raw features, the number of which depends on the instance size. Then we describe how we aggregate these features to obtain a fixed-size vector.

Several features exploit the fact that the graph $G(\ell) = (V(\ell), A(\ell))$ is the state-transition graph of a dynamic program. We assume that at least one state variable (indicating the step in the decision process) is perfectly known, and that the value of this state variable is strictly increasing along a path. As a result, $G(\ell)$ is a layered graph, where each vertex belongs to a layer. We assume that layers are ordered from 0 to the total number of layers. Let $\tau: V(\ell) \to \mathbb{N}$ be the mapping from the vertices to their layer number. Similarly, for any arc $a$, we define $\chi(a) = \tau(v^+(a)) - \tau(v^-(a))$ as the difference between the layers of its endpoints (called layer span of an arc hereafter). We define the span of a constraint $k \in \mathcal{R}$ as the minimum layer interval $[b^k, e^k]$, with $b^k < e^k$, such that for each vertex $u$ in the set $\{v \in V: \tau(v) < b^k \text{ or } \tau(v) > e^k\}$, the resource consumption $r_a^k$ on all outgoing arcs $a \in \mathcal{A}^+(u)$ is zero.

There may be more than one known state variable, depending on the use case. We denote by $\delta^d(v)$ the value of the $d$-th state variable at vertex $v$. The first state variable indicates the layer, i.e., $\tau(v) = \delta^1(v)$. Unlike $\delta^1$, the other state variables $\delta^d$ (for $d > 1$) are not necessarily monotonic along paths.

We organize the features into three categories: instance-based, vertex-based, and path-based. Path-based features focus on the cost and the current resource consumption of the path $p(\ell)$, vertex-based features focus on the last vertex $v^+(p(\ell))$ of the path, and instance-based features involve global information on the problem instance $(G(\ell), \mathcal{R}(\ell))$. We normalize the value of each feature using min-max normalization, mapping it to the interval $[0, 1]$ based on its minimum and maximum possible values. We determine these minimum and maximum values based on the instance $(G(\ell), \mathcal{R}(\ell))$. For each family of features, we report its expression, the minimum and maximum values that it can take, and the number of features in the family.

**Instance-based features** only depend on $(G(\ell), \mathcal{R}(\ell))$. We list them in Table 1. We compute some features directly from $G$: the cost of each arc (GF2), the density of the graph (GF8), and the cost of the arc divided by the number of layers it spans (GF3). The others are related to the Lagrangian relaxation of the problem: Lagrangian costs of the arcs (GF1), Lagrangian multipliers (GF4), initial optimality gap (GF5), ratio of removed arcs (GF6), and ratio of removed vertices (GF7). The value of $UB_{initial}$ is the initial upper bound for the optimal value for the instance obtained from the preprocessing phase, while $LB$ is the value of the Lagrangian dual problem described in Section 3.3.

**Vertex-based features**, summarized in Table 2, are derived from the terminal vertex $v^+(p(\ell))$ of the path. The first feature (VF1) is the layer of the vertex. It reflects the relative depth of the path within the graph, indicating how far along the solution process the path is. A smaller layer suggests that many arcs remain to be appended to reach the sink, while a larger layer implies proximity to completion. The second set of features (VF2(d)) contains the value of each state variable at this
<table><caption>Table 1: Instance-based features</caption><thead><tr><th></th><th>Feature type</th><th>expression</th><th>min</th><th>max</th><th>number</th></tr></thead><tbody><tr><td>GF1</td><td>Lagrangian cost of arcs</td><td>c&#771;<sup>&lambda;*</sup><sub>a</sub></td><td>min{c&#771;<sup>&lambda;*</sup><sub>a</sub>}</td><td>max{c&#771;<sup>&lambda;*</sup><sub>a</sub>}</td><td>|A|</td></tr><tr><td>GF2</td><td>Cost of arcs</td><td>c<sub>a</sub></td><td>min{c<sub>a</sub>}</td><td>max{c<sub>a</sub>}</td><td>|A|</td></tr><tr><td>GF3</td><td>Cost/layer span of arcs</td><td>c<sub>a</sub>/&chi;(a)</td><td>min{c<sub>a</sub>}/&chi;(a)</td><td>max{c<sub>a</sub>}/&chi;(a)</td><td>|A|</td></tr><tr><td>GF4</td><td>Lagrangian multipliers</td><td>&lambda;<sup>*</sup><sub>k</sub></td><td>min{&lambda;<sup>*</sup><sub>k</sub>}, k &isin; R</td><td>max{&lambda;<sup>*</sup><sub>k</sub>}, k &isin; R</td><td>|R|</td></tr><tr><td>GF5</td><td>Initial gap</td><td>|UB<sub>initial</sub> - LB| / (max(UB<sub>initial</sub>, |LB|))</td><td>0</td><td>1</td><td>1</td></tr><tr><td>GF6</td><td>Ratio of removed arcs</td><td>|A<sup>0</sup> - A| / |A|</td><td>0</td><td>1</td><td>1</td></tr><tr><td>GF7</td><td>Ratio of removed vertices</td><td>|V<sup>0</sup> - V| / |V|</td><td>0</td><td>1</td><td>1</td></tr><tr><td>GF8</td><td>Density</td><td>|A| / |V| / (|V| - 1)</td><td>0</td><td>1</td><td>1</td></tr></tbody></table>

vertex, which captures problem-dependent dynamics embedded in the state space. Given a given DPRC(P), the number of state variables *D* remains constant regardless of the instance size.

<table><caption>Table 2: Vertex-based features</caption><thead><tr><th></th><th>Feature type</th><th>expression</th><th>min</th><th>max</th><th>number</th></tr></thead><tbody><tr><td>VF1</td><td>Layer</td><td>&tau;(v)</td><td>0</td><td>&tau;(t)</td><td>1</td></tr><tr><td>VF2(d)</td><td>Value of state variables d > 1</td><td>&delta;<sup>d</sup>(v)</td><td>min{&delta;<sup>d</sup>(u) }<br>u &isin; V</td><td>max{&delta;<sup>d</sup>(u) }<br>u &isin; V</td><td>D - 1</td></tr></tbody></table>

**Path-based features**, listed in Table 3, depend on the path prefix $p(l)$. The first three features are based on quantities previously defined: the average cost of an arc in the path (PF1), the total cost divided by the current layer (PF2), and the Lagrangian lower bound for any path with prefix *p* (PF3). To define the fourth feature type, let $\bar{r}_{min}^k$ and $\bar{r}_{max}^k$ denote the lower and upper bounds on the possible consumption of resource *k* for any feasible path in $G(l)$. The fourth feature type is the total resource consumption along the path (PF4), and the fifth is the resource consumption divided by the current layer (PF5).

The last three features, (PF6), (PF7), and (PF8), are derived from the following quantities: $\varphi^k(p)$, a lower bound on the cost to complete $p(l)$ into an $s-t$ path satisfying constraint *k*; and $\kappa_{min}^k(p)$ and $\kappa_{max}^k(p)$, the minimum and maximum arc costs, respectively, when attempting to satisfy constraint *k*. We provide their detailed computation in Appendix D.

<table><caption>Table 3: Path-based features</caption><thead><tr><th></th><th>Feature type</th><th>expression</th><th>min</th><th>max</th><th>number</th></tr></thead><tbody><tr><td>PF1</td><td>Average arc cost</td><td>c(p)/|p|</td><td>min{c<sub>a</sub>}<br>a &isin; A</td><td>max{c<sub>a</sub>}<br>a &isin; A</td><td>1</td></tr><tr><td>PF2</td><td>Cost divided by layer</td><td>c(p)/&tau;(v<sup>+</sup>(p))</td><td>min{c<sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>max{c<sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>1</td></tr><tr><td>PF3</td><td>Lagrangian cost</td><td>&nu;&lambda;(p)</td><td>LB</td><td>UB</td><td>1</td></tr><tr><td>PF4</td><td>Resource consumption of k</td><td>r<sup>k</sup>(p)</td><td>r&#771;<sup>k</sup><sub>min</sub></td><td>r&#771;<sup>k</sup><sub>max</sub></td><td>|R|</td></tr><tr><td>PF5</td><td>Resource consumption k divided by layer</td><td>r<sup>k</sup>(p)/&tau;(v<sup>+</sup>(p))</td><td>min{r<sup>k</sup><sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>max{r<sup>k</sup><sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>|R|</td></tr><tr><td>PF6</td><td>Completion satisfying constraint k divided by layer span</td><td>&varphi;<sup>k</sup>(p)/(&tau;(t) - &tau;(v<sup>+</sup>(p)))</td><td>min{c<sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>max{c<sub>a</sub>}/&chi;(a)<br>a &isin; A</td><td>|R|</td></tr><tr><td>PF7</td><td>Min arc cost to satisfy constraint k</td><td>&kappa;<sup>k</sup><sub>min</sub>(p)</td><td>min(0, min{c<sub>a</sub>})<br>a &isin; A</td><td>max(0, max{c<sub>a</sub>})<br>a &isin; A</td><td>|R|</td></tr><tr><td>PF8</td><td>Max arc cost to satisfy constraint k</td><td>&kappa;<sup>k</sup><sub>max</sub>(p)</td><td>min(0, min{c<sub>a</sub>})<br>a &isin; A</td><td>max(0, max{c<sub>a</sub>})<br>a &isin; A</td><td>|R|</td></tr></tbody></table>
### 4.2.2 Feature aggregation

Once the features are normalized, we aggregate feature types whose sizes depend on the instance to ensure that the feature vector is of fixed size. We need to decide which original features to gather (clustering) and how to produce a single value for each cluster. An aggregation uses three parameters: the number of clusters, the clustering type, and the aggregator applied to each cluster. Table 4 summarizes the aggregation techniques we used for each feature type.

Let $q$ denote the number of clusters for a given feature type. To simplify the exposition, when only a single cluster is used ($q = 1$), we apply a dummy clustering method referred to as "all," in which all features of the given type belong to a single cluster.

For constraint-related feature types, we aggregate features for each family of resource constraints $f \in \mathcal{F}$ (see Section 3.1 for the definition of $\mathcal{F}$). For each constraint, $q$ partitions are computed once per instance and shared across all paths within that instance. We consider two clustering methods for feature aggregation: the *q-quantile* and the *q-max-disjoint* methods.

The *q-quantile* clustering method first sorts the $m$ constraints in decreasing order based on a specific criterion, and then partitions them into $q$ clusters such that the first $m/q$ constraints are assigned to the first cluster, the next $m/q$ to the second, and so on. The sorting criteria we use are the Lagrangian multipliers (*q-quant-lag*), the span length $e^k - b^k$ (*q-quant-spanL*), and the maximum arc costs with non-zero resource consumption (*q-quant-maxcost*).

The *q-max-disjoint* clustering method aims to minimize the number of constraints whose spans overlap within a cluster. For this purpose, the algorithm initializes $q$ empty clusters, and assigns constraints cluster by cluster. It initializes each cluster with the remaining constraint that has the minimum span end. Then, the algorithm greedily adds the remaining constraint with the smallest span end among those whose span does not overlap with the previously added constraint. This process continues until no such constraint exists. The procedure is then applied to the next cluster, up to the $q$-th cluster. The constraints that remain unassigned after this phase are assigned one by one to clusters in a *round-robin fashion*, cycling back to the first cluster after using the last one.

We apply distinct aggregation techniques to each cluster of features. In total, we use six aggregation functions, tailored to handle different types of features. The first five are classical aggregation functions: maximum, mean value, standard deviation (std), skewness (skew), and kurtosis (kurt). The last two, skew and kurt, capture third and fourth moments of the feature distributions, providing information about asymmetry and tail behavior, respectively. Since their values are unbounded in $\mathbb{R}$, we apply a sigmoid function to map them to the interval (0, 1), ensuring numerical stability during training. We designed the last aggregator, *mean-inspan*, for constraint-related feature types. It sums up only the features associated with active constraints (i.e., constraint $k$ for which the layer $\tau(v(\ell)) \in [b^k, e^k]$), then divides the result by the number of active constraints.

<table><caption>Table 4: Feature aggregation</caption><thead><tr><th>Feature types</th><th>Number</th><th>Aggregator</th><th>Clustering</th></tr></thead><tbody><tr><td>GF1 GF2 GF3</td><td>|A|</td><td>mean, std, skew, kurt</td><td>all</td></tr><tr><td>GF4</td><td>|R|</td><td>mean, std, skew, kurt</td><td>all, q-quant-lag</td></tr><tr><td>PF4 PF5</td><td>|R|</td><td>mean-inspan</td><td>all, q-quant-lag, q-quant-spanL, q-quant-maxcost, q-max-disjoint</td></tr><tr><td>PF6</td><td>|R|</td><td>max</td><td>all</td></tr><tr><td>PF7 PF8</td><td>|R|</td><td>mean-inspan</td><td>all</td></tr></tbody></table>

### 4.2.3 Supervised learning using the mean square error

We trained the FNN using Supervised Learning. We split our instances in the set $\mathcal{I}$ into a set of training instances $\mathcal{I}_{train}$ and a set of validation instances $\mathcal{I}_{val}$. Let $\mathcal{L}$ be the set of (instance-path pair features, score) pairs (see Section 4.1) collected from all instances in $\mathcal{I}$. We denote by $\mathcal{L}_{train} \subseteq \mathcal{L}$
and $\mathcal{L}_{\text{val}} \subseteq \mathcal{L}$ the subsets of (instance-path pair features, score) pairs collected from $\mathcal{I}_{\text{train}}$ and $\mathcal{I}_{\text{val}}$, respectively.

Let $L_{\Upsilon}(\omega)$ be the mean-square error for the NN parameter $\omega$ on the set $\Upsilon \subseteq \mathcal{L}$. We have:

$$L_{\Upsilon}(\omega) = \frac{1}{|\Upsilon|} \sum_{(\phi(\ell), y(p(\ell))) \in \Upsilon} (\tilde{\psi}(\phi(\ell), \omega) - \hat{y}(p(\ell)))^2,$$

where $\hat{y}(p(\ell)) = \frac{y(p(\ell)) - \text{LB}}{\text{UB}_{\text{initial}} - \text{LB}}$ is the normalized value of $y(p(\ell))$, which helps stabilize learning.

The learning problem is $\min_{\omega} L_{\mathcal{L}}(\omega)$. We solved it using a stochastic gradient descent algorithm that updates $\omega$ based on the gradient of $L_{\mathcal{L}_{\text{train}}}$ at $\omega$. To avoid overfitting, the testing phase uses the value of $\omega$ that gives the minimal value of $L_{\mathcal{L}_{\text{val}}}(\omega)$.

# 5 Use cases

This section presents two CO problems that we use to test ML-ABSM: the single machine total weighted tardiness ($1 \parallel \sum_i w_i T_i$) problem and the temporal knapsack problem (TKP). The DPRC formulation of the two problems is presented in E. For each problem, we provide a formal definition of the problem and then describe the process for generating instances. In particular, we describe how we generate non-trivial instances, i.e. instances for which the value of the Lagrangian problem is not equal to the primal bound given by running BSM with $\Psi = ((\nu^{\lambda^*}, \beta^*, 0))$ without rollout and with a small $\beta^*$.

## 5.1 Single machine total weighted tardiness

### Problem definition<fnref n="1" />

The problem $1 \parallel \sum_i w_i T_i$ is a classical one machine scheduling problem (see Pinedo 2022). Let $\mathcal{J}$ be a set $n$ of jobs. Each job $j \in \mathcal{J}$ has a processing time $p_j \in \mathbb{N}_+$, a due date $d_j \in \mathbb{N}_+$, and a weight $w_j \in \mathbb{N}_+$ representing the importance of the job $j$. Suppose that the job $j$ starts at time $t_j$, it is then processed until time $t_j + p_j$, and its tardiness $T_j$ is defined by $T_j = \max(0, t_j + p_j - d_j)$.

### Instance generation

We used a generator from the literature (documented in OR-library). The generator builds instances given a triplet $(n, \text{TF}, \text{RDD})$, where TF and RDD take values in $\{0.2, 0.4, 0.6, 0.8, 1.0\}$. To generate non-trivial instances, we proceeded as follows. For a fixed seed, a given size $n$, and a (TF, RDD) pair, we iteratively generate instances until five non-trivial instances are generated or the maximum number of 100 iterations is reached. In the latter case, we create fewer than five non-trivial instances.

## 5.2 Temporal knapsack

### Problem definition

The temporal knapsack problem (TKP) (Bartlett et al. 2005) is a generalization of the knapsack problem, where each item is associated with a time window during which it can be added to the knapsack. Let $W$ be the knapsack capacity. Let $\mathcal{I} = \{1, \dots, n\}$ be a set of $n$ items. Each item $i$ has a profit $p_i \in \mathbb{R}^+$, a weight $w_i \in \mathbb{N}$, and a time interval $[s_i, f_i)$, where $s_i, f_i \in \mathbb{N}$, $s_i < f_i$. A feasible solution is a subset $\mathcal{J}$ of $\mathcal{I}$ such that for any time unit $t \in \mathbb{N}$, the sum of the weights of the items

<footnote n="1">¹We chose to use the classical notations for this problem, even if they overlap with notations from other sections, since it does not create any ambiguity.</footnote>
in $\mathcal{J}$ whose time interval contains $t$ is less than or equal to $W$. The objective of TKP is to find a feasible solution $\mathcal{J}$ that maximizes the total profit. Early works include approximation algorithms (Calinescu et al. 2002; Chen, Hassin, and Tzur 2002), dynamic programs (Chen, Hassin, and Tzur 2002; Caprara, Furini, and Malaguti 2013), and a Dantzig-Wolfe decomposition (Caprara, Furini, and Malaguti 2013).

## Instance generation

We generated instances following the description of the second class of instances introduced in (Caprara, Furini, and Malaguti 2013), with a slight modification to vary the knapsack capacity. The profit $p_i$ and weight $w_i$ of each item $i$ are independently drawn from a uniform distribution over the interval $[1, 100]$. To generate the time interval $[s_i, f_i)$, the start time $s_i$ is drawn uniformly from the interval $[1, 1000]$, and the duration $f_i - s_i$ is drawn uniformly from $[1, d_{max}]$ with $d_{max} \in \{60, 70, 80, 90, 100\}$. The knapsack capacity $W$ is set to $\alpha \sum_{i=1}^{n} w_i$ with $\alpha \in \{0.01, 0.005, 0.002\}$. We generated instances for a given size $n$ and a given $\{d_{max}, \alpha\}$ parameter pair. Preliminary experimental results suggest that none of the resulting TKP instances were trivial.

# 6 Computational experiments

We now provide numerical evidence that our ML-augmented ABS enables a clear improvement over classical strategies based solely on lower and upper bounds. We first describe the experimental setting, including the learning setup for our ML models, the instance generation parameters, and the evaluation metrics for ABSM, and then present the computational results for the two problems described above. Our experiments focus on the following three research questions. First, what are the best parameters for the ML-augmented ABSM (ML architecture, feature selection, and the proportion of paths to select with ML)? Second, does ML improve the quality of the solution given by ABSM? Third, for a given problem and setting, is the ML model capable of generalizing to larger instances?

## 6.1 Experimental setting

We investigate the effectiveness of ML in heuristics for obtaining primal bounds within an exact solution algorithm. Therefore, for each problem, we consider two settings: (1) when the vector of optimal Lagrangian multipliers $\lambda^*$ is available (i.e. $\lambda = \lambda^*$), and (2) when $\lambda^*$ is not available (i.e., $\lambda = 0$). The former case typically corresponds to the initialization of an iterative algorithm, when one needs an initial upper bound to help the filtering algorithm. The Lagrangian dual bound of a path $p$ is $\nu^{\lambda^*}(p)$. In the second case, we aim to quickly obtain a feasible solution without computing $\lambda^*$. The Lagrangian dual bound of a path $p$ is $\nu^0(p)$.

For DPRC($1 \mid \sum_i w_i T_i$), we solve the Lagrangian dual problem (4) using column generation to compute $\lambda^*$. As a rollout algorithm, we use the greedy algorithm based on the Earliest Due Date (EDD) rule. To generate non-trivial instances for the learning pipeline, we use $n = 40$, four random seeds, and set $\beta' = 20$. This results in a total of $|\mathcal{I}| = 450$ instances, with $|\mathcal{I}_{train}| = 337$ used for training (70%) and $|\mathcal{I}_{val}| = 113$ for validation (30%). In the beam-search run in the pre-processing phase, we set $\beta_0 = 20$. During the sampling and score collection phases, we use $\beta = 100$, $M = 80$, and initialize $\beta_{avg} = \beta - 1$. The value of $\beta_{avg}$ is then decreased by 10 every 10 iterations. The dynamic adjustment of $\beta_{avg}$ over the iterations is motivated by the need to recalibrate the balance between paths sampled under the guidance of $\psi_{avg}$ and those guided by $\psi_{rand}$. The former favors paths that have historically yielded high-quality solutions or exhibit low Lagrangian bounds, whereas the latter encourages exploration of diverse paths. Consequently, decreasing the influence of $\psi_{avg}$ enhances diversification. The enhancement phase requires a beam search evaluation on each sample path, which is computationally heavy; thus, we set $N_{enhance} = 4000$ and $\beta_{enhance} = 10$. We provide
the complete list of features in Table 14. For training, we use an FNN with four hidden layers, each consisting of 200 neurons with ReLU activation functions between consecutive layers. The maximum number of training epochs is 3000, and the batch size is 16384(= 2¹⁴). The initial learning rate is set to 0.05 and is reduced by a factor of 0.1 if the validation loss $L_{\mathcal{L}_{\text{val}}}$ does not improve for 10 consecutive epochs, down to a minimum of 10⁻⁶.

For DPRC(TKP), the Lagrangian dual problem (4) is solved using a fine-tuned volume algorithm (Barahona and Anbil 2000), which proved faster in preliminary experiments than the column generation method for TKP. As a rollout algorithm, we use a greedy algorithm that considers items in increasing order of their start time, selecting the next item only if the knapsack constraint remains satisfied. To generate instances for the learning pipeline, we set $n = 1000$, and generate 40 instances for each $d_{\max}/\alpha$ configuration. This yields a total of $|\mathcal{I}| = 600$ instances, which are randomly split into two subsets: $|\mathcal{I}_{\text{train}}| = 420$ for training (70%) and $|\mathcal{I}_{\text{val}}| = 180$ for validation (70%). In the preprocessing phase, we again set $\beta_0 = 20$. The sampling and score collection phase follows the same settings as for $1||\sum_i w_i T_i$, except that we set $\beta$ to 50. We provide the complete feature list in Table 14. For training, we adopt the same settings as in the $1||\sum_i w_i T_i$ case, with the exception that the FNN uses four hidden layers of 1000 neurons each, with ReLU activation functions between consecutive layers.

The parameters of the testing phase are as follows. For DPRC($1||\sum_i w_i T_i$), we considered two instance sizes: 100 and 200. For each size, we generated 110 instances (220 test instances overall). We used $\beta' = 50$ using one random seed, to check the instance triviality quickly. For DPRC(TKP), we generated 150 test instances for $n = 1000$ and for $n = 2000$, with 10 instances per ($d_{\max}$, $\alpha$) setting, yielding 300 instances in total. In the pre-processing stage of the testing phase, we first executed a BSM with $\Psi = ((\nu^\lambda, \beta_0, 0))$ without rollout, and using small beam sizes ($\beta_0 = 50$ for DPRC($1||\sum_i w_i T_i$) and $\beta_0 = 20$ for DPRC(TKP)) to compute an initial primal bound $\text{UB}_{\text{initial}}$. The initial beam size of ABSM matches $\beta_0$ from the preprocessing phase.

For each use case, we compare the solution quality produced by ABSM against a problem-specific heuristic algorithm, denoted by heur. For $1||\sum_i w_i T_i$, heur corresponds to a fine-tuned *iterated local search* (ILS) (detailed in F). For TKP, heur is *CPLEX* with a time limit of 5 minutes (CPLEX-5min). We evaluate the solution quality using the metric: $\text{gap} = \frac{\text{UB}_{\text{ABSM}} - \text{UB}_{\text{heur}}}{|\text{UB}_{\text{heur}}|}$, where $\text{UB}_{\text{heur}}$ denotes the objective value of the solution returned by heur, and $\text{UB}_{\text{ABSM}}$ is the objective value of the solution obtained by ABSM. The gap can be negative if ABSM returns a better solution than heur.

We implemented the data collection and testing phases in C++ (GCC 12.4.0, glibc 2.39) and trained the neural network using Python 3.12.1 and PyTorch 2.4.0. In the column generation procedure, we used CPLEX 22.1.1. to solve the master problem. We use the same solver as a heuristic for TKP. We conducted the sampling and testing phases on CPU-based machines equipped with 18 cores (Intel Xeon Gold 6240 @ 2.60GHz) and 192GB of RAM. During testing, we ran three tasks in parallel on each machine. We conducted the learning phase on a GPU-based machine equipped with an NVIDIA V100 GPU and 384GB of RAM. To enable fast inference on CPU, we applied dynamic quantization<fnref n="2" /> to the FNN models, reducing the weights from `float32` to `float16`. This approach preserves prediction accuracy while decreasing inference time. We also tried a quantization to `int8`, but it considerably degraded the prediction quality.

In the following subsections 6.2, 6.3, 6.4, we present the numerical results of ML-augmented ABSM on the two problems described above: DPRC($1||\sum_i w_i T_i$) and DPRC(TKP). We report, for each problem, the aggregated results of our two possible settings, with optimal ($\lambda^*$) or null (0) Lagrangian multipliers. Complete data and detailed results are available on gitlab<fnref n="3" />. For clarity, we denote ABSM with ML integration as ML-ABSM, and ABSM relying solely on the Lagrangian lower bound and the rollout upper bound as ABSM. ML scoring is computationally more expensive than traditional scoring due to feature computation and inference overhead. To provide a fair comparison between ML-ABSM and ABSM, we give the same time limit (300 seconds) to all methods.

<footnote n="2">²For background on quantization, see https://www.maaartengrootendorst.com/blog/quantization/</footnote>

<footnote n="3">³https://gitlab.inria.fr/edge/dprcmlbs_results</footnote>
First, we identify the optimal proportion of paths to select with each scoring function (with and without ML), independently of computational cost in Section 6.2. Second, we balance prediction quality and computational efficiency by reducing the feature set and adjusting the ML model size in Section 6.3. Third, by combining lightweight ML models with the optimal path-selection proportions determined in the first part, we analyze how solution quality evolves over time, comparing ML-ABSM to ABSM in Section 6.4.

## 6.2 Tuning the path selection ratios

In this section, we determine the best configuration for the proportion of paths selected by each of the three scoring functions: the Lagrangian dual (LB), the greedy rollout algorithm (Roll), and the ML model. The ratio of paths allocated to each of these functions defines a configuration. We conducted experiments with an ML model trained on the collected data, with feature selection and model size tuning determined from preliminary experiments.

First, we determined the best-performing configuration without ML. We obtained this configuration by testing different proportions of paths selected by the Lagrangian lower bound (rLB) from the set {0, 0.1, ..., 0.9, 1.0}, with the remaining proportion of paths selected by the upper bound given by the rollout algorithm (rRoll), i.e., 1 - rLB. We report the selected feature set in Table 14, and the details of our experiments in G. We also report the best configurations (rLB, rRoll) in Table 5. Unless stated otherwise, we used these configurations in ABSM. We observe that switching from $\lambda^*$ to 0 increases the proportion of paths to be scored by the greedy rollout algorithm, as LB gives by definition worst lower bounds, and therefore provides less guidance toward high-quality paths. The increase is more pronounced for DPRC(1||$\sum_i w_i T_i$) than for DPRC(TKP), as the greedy rollout produces better solutions in the former case.

<table><caption>Table 5: Best configuration for ABSM</caption><thead><tr><th>DPRC(P)</th><th>λ</th><th>rLB</th><th>rRoll</th></tr></thead><tbody><tr><td rowspan="2">1||$\sum_i w_i T_i$</td><td>λ*</td><td>0.80</td><td>0.20</td></tr><tr><td>0</td><td>0.10</td><td>0.90</td></tr><tr><td rowspan="2">TKP</td><td>λ*</td><td>1.00</td><td>-</td></tr><tr><td>0</td><td>0.80</td><td>0.20</td></tr></tbody></table>

Second, we varied the proportion of paths selected by the ML model (rML) from the set {0, 0.25, 0.5, 0.75, 1.0}, and distributed the remaining paths between LB and Roll according to the ratio rLB/rRoll of the previously determined best configuration for each case. We report the selected model size in Table 15. We provide detailed results in H.2 and report the best configurations (rLB, rRoll, rML) in Table 6. In the following, we adopt these configurations for ML-ABSM. We observe that the integration of ML into ABSM is more pronounced for DPRC(1||$\sum_i w_i T_i$) than for DPRC(TKP), suggesting that in the former case the ML model more frequently selects paths that are ignored by the non-ML scoring functions and that lead to better solutions than those retained by the latter.

<table><caption>Table 6: Best configuration for ML-ABSM</caption><thead><tr><th>DPRC(P)</th><th>λ</th><th>rLB</th><th>rRoll</th><th>rML</th></tr></thead><tbody><tr><td rowspan="2">1||$\sum_i w_i T_i$</td><td>λ*</td><td>0.20</td><td>0.05</td><td>0.75</td></tr><tr><td>0</td><td>0.05</td><td>0.45</td><td>0.50</td></tr><tr><td rowspan="2">TKP</td><td>λ*</td><td>0.75</td><td>-</td><td>0.25</td></tr><tr><td>0</td><td>0.60</td><td>0.15</td><td>0.25</td></tr></tbody></table>
## 6.3 Feature selection and model size tuning

This section presents the feature selection and model size tuning processes used to obtain a lightweight ML model. We first conduct a feature analysis for each use case, with the model size fixed as reported in Table 15. Following Parmentier and T'Kindt (2023), we adopted a backward feature elimination method to identify the most relevant features. The process is iterative: at each step, it removes the features ranked as least important by SHAP, retrains the model, and evaluates it by running ML-ABSM on a small validation set of unseen instances with 200 jobs for $1 || \sum_i w_i T_i$ and 2000 items for TKP. It stops when it observes a significant performance degradation of ML-ABSM. We report the detailed results in H.3, and the final feature set in Table 21. For the case DPRC(TKP) with $\lambda = \lambda^*$, the integration of ML into ABSM does not improve the solution quality. The feature set is reduced to a single element, namely the Lagrangian lower bound. This indicates that, in this case, the Lagrangian lower bound (LB) always dominates the ML scoring, in the sense that the best computed solutions correspond to paths whose subpaths from $s$ have all been selected by LB. In other words, ML fails to identify potential high-quality paths that were overlooked by LB. After fixing the feature set, we tuned the model size by varying two parameters: the number of hidden layers and the number of neurons per hidden layer. We tested FNNs with 2, 3, or 4 hidden layers and 50, 100, or 200 neurons per layer, resulting in a total of nine settings. These choices ensure compact ML models with fast inference times. We evaluated each setting by running ML-ABSM on a small validation set of unseen instances. The results are reported in Table 7, where the best-performing model sizes are highlighted in bold. The table also includes the number of paths scored by ML per second. This metric accounts for both feature computation and model inference.

As expected, for a fixed use case, the number of paths scored per second decreases as the model size increases, since the computational load becomes heavier. For a fixed network architecture (with the same number of hidden layers and neurons per layer), the scoring speed varies across use cases due to the differing feature sets. Generally, a larger number of features results in more intensive computations. Notably, for DPRC(TKP) with $\lambda = 0$, feature computation is the most expensive among the four cases, even though this setting does not involve the largest feature set. This is because many features in this case require $O(|R|)$ computations.

In summary, reducing both the feature set and the model size speeds up the scoring of paths by the ML model. In the remainder of these experiments, we use the lightweight ML model to perform the comparisons.

## 6.4 Comparison between ML-ABSM and ABSM

In this section, we assess and compare the evolution of solution quality over time for ABSM with and without ML, specifically examining the impact of ML integration. Our main goal is to compare the objective value obtained by the two variants over time. Average gaps alone do not give enough information to understand the dynamics of the algorithm. Thus, for each tested instance size, we report four figures. The first (leftmost) plot illustrates five curves representing the mean gap(%) of ML-ABSM and ABSM, both under their best configurations and when using only one of the three scoring functions. The curves are shown as a function of computation time, providing an averaged comparison of the relative solution quality over time. To complement it, the second figure depicts the temporal evolution of the distribution of instance-wise relative gap differences between the best solutions obtained so far by ML-ABSM and ABSM, both using their respective best configurations. The mean gap difference is represented by a solid line, with dashed lines, and the shaded region corresponds to the 5th-95th percentile interval. The third figure reports, for each time point, the proportion of instances where ML-ABSM yields superior, inferior, or identical performance relative to the version without ML. Since this representation does not convey the magnitude of performance deviations, the fourth (rightmost) figure presents a per-instance comparison of the relative gaps at the time limit (300 seconds). This enables the quantification of improvements on an instance-by-instance basis. Together, these figures offer a comprehensive overview of the results obtained with
<table><caption>Table 7: Performance of ML-ABSM after 300 seconds using ML models of different sizes</caption><thead><tr><th rowspan="2">DPRC(P)</th><th rowspan="2">λ</th><th rowspan="2">#feat.</th><th rowspan="2">#neu.</th><th rowspan="2">#hid.</th><th colspan="2">Gap(%)</th><th colspan="2">#paths/s</th></tr><tr><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td rowspan="20">TKP</td><td rowspan="6">λ*</td><td rowspan="6">1</td><td>50</td><td>2</td><td><strong>-0.016</strong></td><td>0.026</td><td>10226116</td><td>6284283</td></tr><tr><td>50</td><td>3</td><td><strong>-0.016</strong></td><td>0.026</td><td>3180686</td><td>386822</td></tr><tr><td>50</td><td>4</td><td><strong>-0.016</strong></td><td>0.026</td><td>2536091</td><td>351131</td></tr><tr><td>100</td><td>2</td><td><strong>-0.016</strong></td><td>0.026</td><td>2811311</td><td>389254</td></tr><tr><td>100</td><td>3</td><td><strong>-0.016</strong></td><td>0.026</td><td>2049367</td><td>510336</td></tr><tr><td>100</td><td>4</td><td><strong>-0.016</strong></td><td>0.026</td><td>1144045</td><td>85998</td></tr><tr><td rowspan="4">0</td><td rowspan="4">20</td><td>200</td><td>2</td><td><strong>-0.016</strong></td><td>0.026</td><td>1192133</td><td>83859</td></tr><tr><td>200</td><td>3</td><td><strong>-0.016</strong></td><td>0.026</td><td>756534</td><td>74177</td></tr><tr><td>200</td><td>4</td><td><strong>-0.016</strong></td><td>0.026</td><td>480532</td><td>20738</td></tr><tr><td>50</td><td>2</td><td>5.853</td><td>1.393</td><td>44634</td><td>1249</td></tr><tr><td>50</td><td>3</td><td>6.084</td><td>1.413</td><td>44005</td><td>1266</td></tr><tr><td>50</td><td>4</td><td>6.489</td><td>1.817</td><td>44243</td><td>998</td></tr><tr><td>100</td><td>2</td><td>5.844</td><td>1.557</td><td>42684</td><td>744</td></tr><tr><td rowspan="4">0</td><td rowspan="4">20</td><td>100</td><td>3</td><td>6.014</td><td>1.592</td><td>42643</td><td>554</td></tr><tr><td>100</td><td>4</td><td>6.471</td><td>1.682</td><td>42354</td><td>567</td></tr><tr><td>200</td><td>2</td><td>5.933</td><td>1.406</td><td>41891</td><td>638</td></tr><tr><td>200</td><td>3</td><td><strong>5.610</strong></td><td>1.371</td><td>41098</td><td>889</td></tr><tr><td>200</td><td>4</td><td>6.292</td><td>1.528</td><td>40098</td><td>510</td></tr><tr><td rowspan="10">1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td rowspan="6">λ*</td><td rowspan="6">48</td><td>50</td><td>2</td><td>0.815</td><td>2.348</td><td>73854</td><td>41571</td></tr><tr><td>50</td><td>3</td><td>0.820</td><td>2.593</td><td>72659</td><td>40790</td></tr><tr><td>50</td><td>4</td><td>0.869</td><td>2.290</td><td>71502</td><td>40609</td></tr><tr><td>100</td><td>2</td><td>0.950</td><td>2.538</td><td>71936</td><td>40705</td></tr><tr><td>100</td><td>3</td><td>0.819</td><td>2.276</td><td>70496</td><td>39544</td></tr><tr><td>100</td><td>4</td><td>0.911</td><td>2.375</td><td>68388</td><td>38344</td></tr><tr><td rowspan="4">0</td><td rowspan="4">11</td><td>200</td><td>2</td><td><strong>0.765</strong></td><td>2.291</td><td>68266</td><td>38647</td></tr><tr><td>200</td><td>3</td><td>0.896</td><td>2.509</td><td>64693</td><td>36114</td></tr><tr><td>200</td><td>4</td><td>1.150</td><td>2.547</td><td>61329</td><td>34287</td></tr><tr><td>50</td><td>2</td><td>6.090</td><td>5.935</td><td>262097</td><td>71289</td></tr><tr><td rowspan="6">0</td><td rowspan="6">11</td><td><strong>50</strong></td><td><strong>3</strong></td><td><strong>5.859</strong></td><td>5.869</td><td>256503</td><td>71144</td></tr><tr><td>50</td><td>4</td><td>5.979</td><td>5.978</td><td>250070</td><td>67692</td></tr><tr><td>100</td><td>2</td><td>5.986</td><td>5.926</td><td>255305</td><td>69806</td></tr><tr><td>100</td><td>3</td><td>5.864</td><td>5.873</td><td>234807</td><td>64279</td></tr><tr><td>100</td><td>4</td><td>5.927</td><td>5.877</td><td>217720</td><td>57585</td></tr><tr><td>200</td><td>2</td><td>6.476</td><td>6.099</td><td>224524</td><td>61223</td></tr><tr><td>200</td><td>3</td><td>5.882</td><td>5.886</td><td>195161</td><td>52270</td></tr><tr><td>200</td><td>4</td><td>6.203</td><td>6.286</td><td>173972</td><td>46622</td></tr></tbody></table>

and without ML and facilitate a detailed comparison.

To simplify the notation, we denote by LB-ABS, Roll-ABS, and ML-ABS, the ABS variants using solely the Lagrangian lower bound, the upper bound from the rollout algorithm, and the ML scoring function, respectively.

Figure 2 summarizes our experiments for DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>) with λ = λ*. For four out of the five variants, the mean gap decreases over time, indicating continuous improvement in solution quality, except for the rollout-only version of ABSM. We first observe that combining multiple scoring functions yields better average performance than relying on a single scoring function. For both instance sizes, the mean gap of ML-ABSM is consistently lower than that of ABSM. The gap difference between the two methods is more significant for n = 200. The mean gap difference is consistently below 0, indicating that ML-ABSM outperforms ABSM on average at any time. The gap difference at the 95th percentile slightly exceeds 0 for n = 200, suggesting that instances where ML-ABSM performs worse experience only minor degradation. In contrast, the gap difference at the 5th percentile is significantly below 0 for both n = 100 and n = 200, demonstrating substantial improvement on instances where ML-ABSM performs better. For numerous instances, the values of the solutions of ML-ABSM and ABSM remain identical over time. This is unsurprising, as the strong Lagrangian lower bound (LB) associated with λ*, combined with the greedy rollout algorithm (Roll), already ensures the selection of high-quality paths. For the remaining instances, ML-ABSM provides substantial improvements over ABSM for a significant proportion of them, while exhibiting only a slight difference in solution quality for the few instances in which it decreases the performance
<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison of ABSM between with and without ML on DPRC($1/\|\sum_i w_i T_i\|$) with $\lambda = \lambda^*$</figcaption></figure>

of beam search. This trend becomes more pronounced as $n$ increases. Overall, this highlights that ML-ABSM is capable of improving ABSM results, even when the Lagrangian lower bound already provides strong guidance toward high-quality solutions. One can observe that the two methods are competitive with an ad-hoc local-search-based heuristic dedicated to the problem. While the dual bound (LB) already serves as an excellent path selection criterion, an ML-learned estimation of an enhanced rollout identifies overlooked paths that lead to high-quality solutions.

Figure 3 reports our experiments for DPRC($1/\|\sum_i w_i T_i\|$) with $\lambda = 0$. Note that in this case, the quality of the Lagrangian lower bound is significantly worse than when $\lambda = \lambda^*$. Again, the mean gap reduces over time for the five variants, demonstrating improvements in solution quality with increasing computation time and beam width. On average, ML-ABSM reaches solutions in a few seconds that are at least as good as those obtained by ABSM after 300 seconds. The superior average performance of ML-ABSM over ABSM is evidenced by the consistently negative mean gap difference. The gap difference at the 95th percentile is slightly greater than 0, indicating that for instances where ML-ABSM performs worse than ABSM, the degradation is minor. Conversely, the 5th percentile gap difference is well below 0, indicating substantial improvements in instances where ML-ABSM performs better. Notably, the third figure indicates that ML-ABSM is better than ABSM in roughly 40% of the instances, but is worse in nearly half of them for $n = 200$. The fourth figure illustrates this seemingly counterintuitive result: only very small decreases in performance occur, while the improvements are significantly more substantial. This behavior can be explained by the guidance provided by the combination of Lagrangian dual bound and greedy rollout. While ABSM results indicate that there is room for improvement due to the inaccuracy of the Lagrangian dual bound as $n$ increases, the ML scoring function can offer weak guidance in some instances because it partially relies on this bound as one of its input features. Even with substantial improvements relative to ABSM, the solutions remain distant from those achieved with $\lambda^*$, emphasizing the importance of the dual bound as a feature. Both the proportion of instances showing improvement and the magnitude of these improvements are slightly higher for the 100-job instances than for the 200-job instances, suggesting that the ML model may generalize less effectively to larger instances of this problem.

In Figure 4, we report our results for DPRC(TKP) with $\lambda = \lambda^*$. The mean-gap curves of variants that include the Lagrangian lower bound as a scoring function almost perfectly overlap.
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of ABSM between with and without ML on DPRC($1/||\sum_i w_i T_i||$) with $\lambda = 0$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of ABSM between with and without ML on DPRC(TKP) with $\lambda = \lambda^*$</figcaption></figure>
Within the first 50 seconds, the solution quality improves dramatically, reaching a very small mean gap. The mean gap difference remains close to 0 throughout the runtime. While the solutions differ for most instances initially, the difference diminishes over time. At the end of 300 seconds, ML-ABSM produces the same solutions as ABSM for all instances with $n = 1000$ and for about 80% of instances with $n = 2000$. The observed behavior is not surprising, as the Lagrangian cost is the only feature used by the ML method. This is a typical case where the added value of ML is negligible, although it does not hurt the overall performance of the method.

<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of ABSM between with and without ML on DPRC(TKP) with $\lambda = 0$</figcaption></figure>

Figure 5 reports results for DPRC(TKP) with $\lambda = 0$. For the five variants, the mean gap decreases over time. On average, ML-ABSM outperforms ABSM. The 5th–95th percentile range lies largely below 0, indicating that, over time, the improvements in solution quality achieved by ML-ABSM for certain instances relative to ABSM consistently exceed any decreases observed for other instances. ML-ABSM yields worse solutions than ABSM for a few instances, with only a minor decrease in solution quality for most. Improvements are observed in roughly two-thirds of the instances, typically larger than the deterioration in solution quality observed for the others. The results in this problem are superior on instances of larger size than those used for training. This is a case where integrating ML improves the method’s performance, although it can occasionally degrade it in specific instances.

To conclude this comparison, we also report the final gaps of the methods in Tables 8 and 9. We report six configurations: pure Roll, pure LB, pure ML, the best configuration without ML, the best configuration with ML, and random scoring. In the random scoring configuration, the best configuration with ML is retained, and the ML scoring is replaced by a random function that assigns scores drawn independently from the uniform distribution $U(0, 1)$. The random test is performed with five independent runs, each using a different random seed for every instance. Comparing our method with the random score helps assess whether the performance of ML is due to random diversification of the solutions, or if the score truly captures the structure of the instance. For $1 \le \sum_i w_i T_i \le 2$, ML-ABSM provides the best gaps on average, followed by ABSM. Note that replacing ML with a random choice leads to significantly weaker results. This is also true for TKP with $\lambda = 0$. On the contrary, for TKP with $\lambda = \lambda^*$, the best version of the method is LB-ABS. In this case, using random paths does not improve the results, although it is slightly better than ML-ABSM.
for $n = 200$. These results suggest that random diversification alone is insufficient to improve the results. They confirm that ML captures information that is useful for path selection.

<table><caption>Table 8: Performance of ABSM variants at the time limit (300 seconds) with lightweight ML models on DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>)</caption><thead><tr><th rowspan="2">DPRC(P)</th><th rowspan="2">λ</th><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th rowspan="2">rRnd</th><th colspan="2">n = 100<br>Gap (%)</th><th colspan="2">n = 200<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="10">1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td rowspan="5">λ*</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>1.211</td><td>2.822</td><td>1.585</td><td>2.615</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>-</td><td>0.420</td><td>1.837</td><td>0.952</td><td>2.028</td></tr><tr><td>-</td><td>-</td><td>1.00</td><td>-</td><td>0.289</td><td>1.687</td><td>0.923</td><td>1.797</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td>-</td><td>0.264</td><td>1.465</td><td>0.747</td><td>1.576</td></tr><tr><td>0.20</td><td>0.05</td><td>0.75</td><td>-</td><td><b>0.177</b></td><td>1.418</td><td><b>0.544</b></td><td>1.377</td></tr><tr><td>0.20</td><td>0.05</td><td>-</td><td>0.75</td><td>0.264</td><td>1.460</td><td>0.743</td><td>1.572</td></tr><tr><td rowspan="5">0</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>-</td><td>5.334</td><td>5.875</td><td>8.558</td><td>7.283</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>-</td><td>-</td><td>13.125</td><td>11.274</td><td>26.915</td><td>20.622</td></tr><tr><td>-</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>4.536</td><td>6.205</td><td>11.238</td><td>13.627</td></tr><tr><td>0.10</td><td>0.90</td><td>-</td><td>-</td><td>-</td><td>5.003</td><td>5.003</td><td>7.734</td><td>6.481</td></tr><tr><td>0.05</td><td>0.45</td><td>0.50</td><td>-</td><td><b>3.323</b></td><td>4.437</td><td><b>6.118</b></td><td>7.100</td><td></td></tr><tr><td>0.05</td><td>0.45</td><td>-</td><td>0.50</td><td>0.50</td><td>5.003</td><td>4.985</td><td>7.702</td><td>6.443</td></tr></tbody></table>

<table><caption>Table 9: Performance of ABSM variants at time limit (300 seconds) with lightweight ML models on DPRC(TKP)</caption><thead><tr><th rowspan="2">DPRC(P)</th><th rowspan="2">λ</th><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th rowspan="2">rRnd</th><th colspan="2">n = 1000<br>Gap (%)</th><th colspan="2">n = 2000<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="10">TKP</td><td rowspan="5">λ*</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>0.172</td><td>0.191</td><td>0.345</td><td>0.173</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>-</td><td><b>-0.001</b></td><td>0.005</td><td><b>-0.011</b></td><td>0.022</td></tr><tr><td>-</td><td>-</td><td>1.00</td><td>-</td><td><b>-0.001</b></td><td>0.005</td><td>-0.010</td><td>0.023</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>-</td><td><b>-0.001</b></td><td>0.005</td><td><b>-0.011</b></td><td>0.022</td></tr><tr><td>0.75</td><td>-</td><td>0.25</td><td>-</td><td><b>-0.001</b></td><td>0.005</td><td>-0.009</td><td>0.024</td></tr><tr><td>0.75</td><td>-</td><td>-</td><td>0.25</td><td><b>-0.001</b></td><td>0.005</td><td><b>-0.011</b></td><td>0.022</td></tr><tr><td rowspan="5">0</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>-</td><td>6.221</td><td>4.028</td><td>15.519</td><td>3.051</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.984</td><td>0.916</td><td>6.559</td><td>1.481</td></tr><tr><td>-</td><td>-</td><td>1.00</td><td>-</td><td>-</td><td>2.060</td><td>1.598</td><td>10.498</td><td>1.760</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td>-</td><td>-</td><td>0.941</td><td>0.881</td><td>6.283</td><td>1.451</td></tr><tr><td>0.60</td><td>0.15</td><td>0.25</td><td>-</td><td><b>0.917</b></td><td>0.840</td><td><b>5.824</b></td><td>1.520</td><td></td></tr><tr><td>0.60</td><td>0.15</td><td>-</td><td>0.25</td><td>0.941</td><td>0.879</td><td>6.278</td><td>1.454</td><td></td></tr></tbody></table>

## 6.5 Summary of the experiments

We observed significant improvements with the integration of an ML scoring function into ABSM for DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>), both with and without optimized Lagrangian multipliers. The LB scoring function, combined with the greedy rollout algorithm, is insufficient to find high-quality solutions. In contrast, the results for DPRC(TKP) are more nuanced. For the case where λ = 0, performance improves on average, with notable gains for a substantial proportion of instances; however, in a few cases, solution quality deteriorates, but in a rather limited way. When λ = λ*, the strong performance of the Lagrangian lower bound leaves little room for ML to complement path selection by identifying overlooked good paths. Overall, despite the computational overhead of the ML scoring function, integrating ML into ABSM consistently improves solution quality on average compared to ABSM without ML. Learning the values obtained by an enhanced rollout algorithm provides additional guidance that almost always benefits the method in our experiments. When ML-ABSM produces solutions of lower quality than ABSM, the decrease is generally minor, indicating that using an ML scoring function is not detrimental to the method and offers benefits that greatly outweigh any limitations. In general, ML-ABSM achieves better solutions earlier and maintains this
advantage throughout the algorithm. Improvements of similar magnitude are observed as instance
size grows, demonstrating the generalizability of the ML model. Altogether, this supports the
systematic integration of ML into ABSM.

7 Conclusion

We proposed a generic Machine Learning (ML)-augmented Beam Search (BS) for solving Combi-
natorial Optimization problems that admit a resource-constrained shortest path reformulation. The
ML model is used to define a scoring function that complements the guidance from lower and upper
bounds during path selection.

We developed a complete learning pipeline, including data collection, supervised training, and
evaluation within anytime BS (ABS). Our approach relies on generic features extracted directly
from the RCSPP-DAG reformulation of CO problems. To ensure compatibility across instances of
varying sizes, we propose a systematic method for aggregating features and providing a fixed-size
input to a Feed-Forward Neural Network.

The integration of ML within ABS with multiple scoring functions has been tested on two problems, $1 \parallel \sum_i w_i T_i$ and TKP, under both settings: with and without access to optimal Lagrangian multipliers. In all cases except for TKP with access to optimal Lagrangian multipliers, integrating an ML-based scoring function into ABSM improves solution quality. This improvement holds for test instances larger than those used during training. Our tests show that combining multiple scoring functions reduces the risk of misguidance that can occur when relying on a single scoring function. In addition, they indicate that an ML scoring function can be integrated into the method in a systematic manner.

Our feature engineering highlights how expert knowledge can simplify the learning problem. By encoding the resource consumption and cost information into compact, aggregated features, the neural network focuses on learning meaningful patterns rather than interpreting the full graph structure. However, this also shifts the representational burden onto the designer: if key information is excluded during feature construction, the ML model may fail to recognize important patterns. Future work could explore replacing manual aggregation with learned embeddings (e.g., potentially using architectures designed to process resource vectors or local subgraph neighborhoods), thereby striking a better balance between expressiveness and complexity.

Lastly, our use of offline supervised learning means the model’s generalization is contingent on the
similarity between training and test distributions. Training on smaller instances is computationally
efficient, but only effective if those instances exhibit the structural patterns that generalize to larger
ones. In scenarios where instance characteristics change significantly with size, it may be necessary
to include larger training instances or adopt fine-tuning strategies.

Acknowledgements

This work was supported by the ANRT (Association nationale de la recherche et de la technologie) with CIFRE contract 2022/1682. The authors would like to express their sincere gratitude to Najib Errami for his support in the development and optimization of the code for this project. His insights and technical assistance significantly contributed to the implementation of the computational framework. The authors also thank Hola Adrakey for his valuable comments and feedback on this paper. Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see http://www.plafrim.fr).
References

Akeb, H., M. Hifi, and R. M'Hallah (2009). "A Beam Search Algorithm for the Circular Packing Problem". In: *Computers & Operations Research* 36.5, pp. 1513–1528.

Avella, P., M. Boccia, and A. Sforza (2002). "A Penalty Function Heuristic for the Resource Constrained Shortest Path Problem". In: *European Journal of Operational Research* 142.2, pp. 221–230.

Barahona, F. and R. Anbil (2000). "The Volume Algorithm: Producing Primal Solutions with a Subgradient Method". In: *Mathematical Programming* 87.3, pp. 385–399.

Bartlett, M., A. M. Frisch, Y. Hamadi, I. Miguel, S. A. Tarim, and C. Unsworth (2005). "The Temporal Knapsack Problem and Its Solution". In: *Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems*. Ed. by R. Barták and M. Milano. Springer, pp. 34–48.

Bengio, Y., A. Lodi, and A. Prouvost (2021). "Machine Learning for Combinatorial Optimization: A Methodological Tour d'horizon". In: *European Journal of Operational Research* 290.2, pp. 405–421.

Bennell, J. A. and X. Song (2010). "A beam search implementation for the irregular shape packing problem". In: *Journal of Heuristics* 16, pp. 167–188.

Blum, C. and M. J. Blesa (2007). "Probabilistic Beam Search for the Longest Common Subsequence Problem". In: *Engineering Stochastic Local Search Algorithms. Designing, Implementing and Analyzing Effective Heuristics*. Ed. by T. Stützle, M. Birattari, and H. H. Hoos. Lecture Notes in Computer Science. Springer, pp. 150–161.

Blum, C. and C. Miralles (2011). "On Solving the Assembly Line Worker Assignment and Balancing Problem via Beam Search". In: *Computers & Operations Research*. Project Management and Scheduling 38.1, pp. 328–339.

Boland, N., J. Dethridge, and I. Dumitrescu (2006). "Accelerated Label Setting Algorithms for the Elementary Resource Constrained Shortest Path Problem". In: *Operations Research Letters* 34.1, pp. 58–68.

Cabrera, N., A. L. Medaglia, L. Lozano, and D. Duque (2020). "An Exact Bidirectional Pulse Algorithm for the Constrained Shortest Path". In: *Networks* 76.2, pp. 128–146.

Calinescu, G., A. Chakrabarti, H. Karloff, and Y. Rabani (2002). "Improved Approximation Algorithms for Resource Allocation". In: *Integer Programming and Combinatorial Optimization*. Ed. by W. J. Cook and A. S. Schulz. Springer, pp. 401–414.

Caprara, A., F. Furini, and E. Malaguti (2013). "Uncommon Dantzig-Wolfe Reformulation for the Temporal Knapsack Problem". In: *INFORMS Journal on Computing* 25.3, pp. 560–571.

Chen, B., R. Hassin, and M. Tzur (2002). "Allocation of Bandwidth and Storage". In: *IIE Transactions* 34.5, pp. 501–507.

Choo, J., Y.-D. Kwon, J. Kim, J. Jae, A. Hottung, K. Tierney, and Y. Gwon (2022). "Simulation-Guided Beam Search for Neural Combinatorial Optimization". In: *Advances in Neural Information Processing Systems* 35, pp. 8760–8772.

Christofides, N., A. Mingozzi, and P. Toth (1981). "State-Space Relaxation Procedures for the Computation of Bounds to Routing Problems". In: *Networks* 11.2, pp. 145–164.

Clautiaux, F., B. Detienne, and G. Guillot (2021). "An iterative dynamic programming approach for the temporal knapsack problem". In: *European Journal of Operational Research* 293.2, pp. 442–456.

Costa, L., C. Contardo, and G. Desaulniers (2019). "Exact Branch-Price-and-Cut Algorithms for Vehicle Routing". In: *Transportation Science* 53.4, pp. 946–985.

Dumitrescu, I. and N. Boland (2003). "Improved Preprocessing, Labeling and Scaling Algorithms for the Weight-Constrained Shortest Path Problem". In: *Networks* 42.3, pp. 135–153.
Ettrich, R., M. Huber, and G. R. Raidl (2023). "A Policy-Based Learning Beam Search for Combinatorial Optimization". In: *Evolutionary Computation in Combinatorial Optimization*. Ed. by L. Pérez Cáceres and T. Stützle. Springer Nature Switzerland, pp. 130–145.

Ghirardi, M. and C. N. Potts (2005). "Makespan Minimization for Scheduling Unrelated Parallel Machines: A Recovering Beam Search Approach". In: *European Journal of Operational Research*. Project Management and Scheduling 165.2, pp. 457–467.

Hart, P. E., N. J. Nilsson, and B. Raphael (1968). "A Formal Basis for the Heuristic Determination of Minimum Cost Paths". In: *IEEE Transactions on Systems Science and Cybernetics* 4.2, pp. 100–107.

Harvey, W. D. and M. L. Ginsberg (1995). "Limited Discrepancy Search". In: *Proceedings of the 14th International Joint Conference on Artificial Intelligence - Volume 1*. IJCAI'95. Morgan Kaufmann Publishers Inc., pp. 607–613.

Himmich, I., H. Ben Amor, I. El Hallaoui, and F. Soumis (2020). "A Primal Adjacency-Based Algorithm for the Shortest Path Problem with Resource Constraints". In: *Transportation Science* 54.5, pp. 1153–1169.

Himmich, I., I. El Hallaoui, and F. Soumis (2024). "A Multiphase Dynamic Programming Algorithm for the Shortest Path Problem with Resource Constraints". In: *European Journal of Operational Research* 315.2, pp. 470–483.

Horváth, M. and T. Kis (2016). "Solving Resource Constrained Shortest Path Problems with LP-based Methods". In: *Computers & Operations Research* 73, pp. 150–164.

Huber, M. and G. R. Raidl (2022). "Learning Beam Search: Utilizing Machine Learning to Guide Beam Search for Solving Combinatorial Optimization Problems". In: *Machine Learning, Optimization, and Data Science*. Ed. by G. Nicosia, V. Ojha, E. La Malfa, G. La Malfa, G. Jansen, P. M. Pardalos, G. Giuffrida, and R. Umeton. Lecture Notes in Computer Science. Springer International Publishing, pp. 283–298.

Ibaraki, T. (1987). "Chapter 7 Successive Sublimation Methods for Dynamic Programming Computation". In: *Annals of Operations Research* 11.1, pp. 397–439.

Kool, W., H. van Hoof, and M. Welling (2019). "Attention, Learn to Solve Routing Problems!" In: *International Conference on Learning Representations*, pp. 1–25.

Kramer, A., E. Lalla-Ruiz, M. Iori, and S. Voß (2019). "Novel Formulations and Modeling Enhancements for the Dynamic Berth Allocation Problem". In: *European Journal of Operational Research* 278.1, pp. 170–185.

Kuroiwa, R. and J. C. Beck (2023). "Solving Domain-Independent Dynamic Programming Problems with Anytime Heuristic Search". In: *Proceedings of the International Conference on Automated Planning and Scheduling* 33, pp. 245–253.

Kwon, Y.-D., J. Choo, I. Yoon, M. Park, D. Park, and Y. Gwon (2021). "Matrix Encoding Networks for Neural Combinatorial Optimization". In: *Proceedings of the 35th International Conference on Neural Information Processing Systems*. NIPS '21. Curran Associates Inc., pp. 5138–5149.

Labassi, A. G., D. Chételat, and A. Lodi (2022). "Learning to Compare Nodes in Branch and Bound with Graph Neural Networks". In: *Proceedings of the 36th International Conference on Neural Information Processing Systems*. NIPS '22. Curran Associates Inc., pp. 32000–32010.

Lee, J., Y. Lee, J. Kim, A. Kosiorek, S. Choi, and Y. W. Teh (2019). "Set Transformer: A Framework for Attention-based Permutation-Invariant Neural Networks". In: *Proceedings of the 36th International Conference on Machine Learning*. PMLR, pp. 3744–3753.

Libralesso, L., P. A. Focke, A. Secardin, and V. Jost (2022). "Iterative beam search algorithms for the permutation flowshop". In: *European Journal of Operational Research* 301.1, pp. 217–234.

Lowerre, B. T. (1976). "The Harpy Speech Recognition System". PhD thesis. Carnegie-Mellon University.

Lozano, L. and A. L. Medaglia (2013). "On an Exact Method for the Constrained Shortest Path Problem". In: *Computers & Operations Research* 40.1, pp. 378–384.
Mazyavkina, N., S. Sviridov, S. Ivanov, and E. Burnaev (2021). "Reinforcement Learning for Combinatorial Optimization: A Survey". In: *Computers & Operations Research* 134, p. 105400.

Narita, M., R. Kuroiwa, and J. C. Beck (2025). "Reinforcement Learning-based Heuristics to Guide Domain-Independent Dynamic Programming". In: *International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research*. Ed. by G. Tack. Springer. Springer Nature Switzerland, pp. 137-154.

Negrinho, R., M. R. Gormley, and G. J. Gordon (2018). "Learning Beam Search Policies via Imitation Learning". In: *Proceedings of the 32nd International Conference on Neural Information Processing Systems*. NIPS'18. Curran Associates Inc., pp. 10675-10684.

Parmentier, A. and V. T'Kindt (2023). "Structured Learning Based Heuristics to Solve the Single Machine Scheduling Problem with Release Times and Sum of Completion Times". In: *European Journal of Operational Research* 305.3, pp. 1032-1041.

Pessoa, A., E. Uchoa, M. P. de Aragão, and R. Rodrigues (2010). "Exact Algorithm over an Arc-Time-Indexed Formulation for Parallel Machine Scheduling Problems". In: *Mathematical Programming Computation* 2.3, pp. 259-290.

Pinedo, M. L. (2022). *Scheduling: Theory, Algorithms, and Systems*. 6th. Springer Cham.

Pugliese, L. D. P. and F. Guerriero (2013). "A Survey of Resource Constrained Shortest Path Problems: Exact Solution Approaches". In: *Networks* 62.3, pp. 183-200.

Thomas, B. W., T. Calogiuri, and M. Hewitt (2019). "An Exact Bidirectional A* Approach for Solving Resource-constrained Shortest Path Problems". In: *Networks* 73.2, pp. 187-205.

Vadlamudi, S. G., S. Aine, and P. P. Chakrabarti (2013). "Anytime Pack Heuristic Search". In: *Pattern Recognition and Machine Intelligence*. Ed. by P. Maji, A. Ghosh, M. N. Murty, K. Ghosh, and S. K. Pal. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, pp. 628-634.

Vadlamudi, S. G., P. Gaurav, S. Aine, and P. P. Chakrabarti (2012). "Anytime Column Search". In: *AI 2012: Advances in Artificial Intelligence*. Ed. by M. Thielscher and D. Zhang. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, pp. 254-265.

Wu, L., P. Cui, J. Pei, and L. Zhao, eds. (2022). *Graph Neural Networks: Foundations, Frontiers, and Applications*. Springer Nature Singapore.

Yilmaz, K. and N. Yorke-Smith (2021). "A Study of Learning Search Approximation in Mixed Integer Branch and Bound: Node Selection in SCIP". In: *AI 2.2 (2)*, pp. 150-178.

Zaheer, M., S. Kottur, S. Ravanbhakhsh, B. Póczos, R. Salakhutdinov, and A. J. Smola (2017). "Deep Sets". In: *Proceedings of the 31st International Conference on Neural Information Processing Systems*. NIPS'17. Curran Associates Inc., pp. 3394-3404.

Zhang, W. (1998). "Complete Anytime Beam Search". In: *Proceedings of the Fifteenth National/Tenth Conference on Artificial Intelligence/Innovative Applications of Artificial Intelligence*. AAAI '98/IAAI '98. USA: American Association for Artificial Intelligence, pp. 425-430.
# A Mathematical models

Here, we present a mixed integer linear programming (MILP) formulation of DPRC(P) and we subsequently apply Lagrangian relaxation to the resource constraints.

## A.1 MILP formulation of DPRC(P)

The MILP formulation of DPRC(P) requires, for each arc $a \in A$, a binary variable, denoted by $x_a$. The variable $x_a$ is equal to 1 if arc $a \in A$ is selected and 0 otherwise. The MILP formulation is as follows:

$$ \min \sum_{a \in A} c_a x_a \qquad (3a) $$

$$ \text{s.t.} \quad \sum_{a \in A^-(v)} x_a - \sum_{a' \in A^+(v)} x_{a'} = \begin{cases} -1, & \text{if } v = s, \\ 1, & \text{if } v = t, \\ 0, & \text{otherwise,} \end{cases} \quad \forall v \in V, \qquad (3b) $$

$$ \sum_{a \in A} r_a^k x_a = R^k, \qquad \forall k \in \mathbb{R}^+, \qquad (3c) $$

$$ \sum_{a \in A} r_a^k x_a \leq R^k, \qquad \forall k \in \mathbb{R}^+, \qquad (3d) $$

$$ x_a \in \{0, 1\}, \qquad \forall a \in A. \qquad (3e) $$

The objective (3a) of RCSPP is to minimize the total cost of arcs. The constraint (3b) is the flow conservation constraint. The constraints (3c)-(3d) are resource constraints.

## A.2 Lagrangian relaxation

We now precisely describe how we perform a Lagrangian relaxation of the resource constraints. Let $\lambda$ be a vector of size $|\mathcal{R}|$, in which each element $\lambda^k \in \mathbb{R}$, $\forall k \in \mathcal{R}^+$, $\lambda^k \in \mathbb{R}^+$, $\forall k \in \mathcal{R}^\le$ represents the Lagrangian multiplier associated with the k-th resource constraint. The Lagrangian problem is

$$ \max_{\lambda} \mathrm{LR}(\lambda), \qquad (4) $$

where

$$ \mathrm{LR}(\lambda) = \min \sum_{a \in A} c_a x_a + \sum_{k \in \mathcal{R}} \lambda^k \left( \sum_{a \in A} r_a^k x_a - R^k \right) \qquad (5a) $$

$$ \text{s.t.} \quad \sum_{a \in A | v^+(a) = v} x_a - \sum_{a' \in A | v^-(a') = v} x_{a'} = \begin{cases} -1, & \text{if } v = s, \\ 1, & \text{if } v = t, \\ 0, & \text{otherwise,} \end{cases} \quad \forall v \in V, \qquad (5b) $$

$$ x_a \in \{0, 1\}, \quad \forall a \in A \qquad (5c) $$

The objective function can be reformulated as $\sum_{k \in \mathcal{R}} \lambda^k R^k + \sum_{a \in A} \tilde{c}_a^\lambda x_a$.

# B Details on BSM

Algorithm 1 uses Algorithm 3 and Algorithm 4 to check the feasibility of the path and remove dominated paths, respectively.
## B.1 Feasibility

To check whether a path $p$ is feasible, we implemented Algorithm 3. This algorithm uses quantities that are computed during preprocessing. Let $r_{\min}^k(v)$ (resp. $r_{\max}^k(v)$) be the minimum (resp. maximum) of the resource consumption $k \in \mathcal{R}$ of $v$-t paths in $\mathcal{G}$. The computation of $r_{\min}^k(v)$ (resp. $r_{\max}^k(v)$) is achieved by determining the shortest (resp. longest) path from $v$ to $t$ on the graph, with the cost of $a \in \mathcal{A}$ assigned to $r_a^k$. Given an instance $(\mathcal{G}, \mathcal{R})$ and a path $p$, the algorithm checks the necessary conditions for the existence of an extension of $p$ into an $s-t$ path satisfying the resource consumption.

<caption>**Algorithm 3:** IsFeasible($\mathcal{G}, \mathcal{R}, p$)</caption>

**Input** : an instance $(\mathcal{G}, \mathcal{R})$ of the problem, a path $p$
**Output** : True if the path may be part of a feasible $s-t$ path satisfying the resource constraints $\mathcal{R}$, False otherwise
1 **for** $k \in \mathcal{R}^\le$ **do**
2 &nbsp;&nbsp;&nbsp;**if** $r^k(p) + r_{\min}^k(v(p)) > R^k$ **then**
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** False;
4 **for** $k \in \mathcal{R}^=$ **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $r^k(p) + r_{\min}^k(v(p)) > R^k$ or $r^k(p) + r_{\max}^k(v(p)) < R^k$ **then**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** False;
7 **return** True;

## B.2 Dominance

Algorithm 1 uses Algorithm 4 to remove dominated paths. This algorithm takes as input an instance $(\mathcal{G}, \mathcal{R})$ and a set $\mathcal{B}$ of paths, and returns a set $\mathcal{B}'$ of non-dominated paths. A path $p \in \mathcal{B}$ is included in $\mathcal{B}'$ if it is not dominated by any other path $p' \in \mathcal{B}$.

We use sufficient conditions for dominance: if $v^+(p') = v^+(p)$, $c(p') \le c(p)$, $r^k(p') \le r^k(p)$ for all $k \in \mathcal{R}^\le$, and $r^k(p') = r^k(p)$ for all $k \in \mathcal{R}^=$, then $p$ is dominated by $p'$.

<caption>**Algorithm 4:** RemoveDominatedPaths($\mathcal{G}, \mathcal{R}, \mathcal{B}$)</caption>

**Input** : an instance $(\mathcal{G}, \mathcal{R})$ of the problem, a set $\mathcal{B}$ of paths
**Output** : a set of non-dominated paths
1 $\mathcal{B}' \leftarrow \emptyset$;
2 **for** $p \in \mathcal{B}$ **do**
3 &nbsp;&nbsp;&nbsp;dominated $\leftarrow$ True;
4 &nbsp;&nbsp;&nbsp;**for** $p' \in \mathcal{B} \setminus \{p\}$ : ($v^+(p) = v^+(p')$) $\land$ ($c(p') \le c(p)$) **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $k \in \mathcal{R}^\le$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $r^k(p') > r^k(p)$ **then** dominated $\leftarrow$ False;
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $k \in \mathcal{R}^=$ **do**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $r^k(p') \neq r^k(p)$ **then** dominated $\leftarrow$ False;
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** not dominated **then** $\mathcal{B}' \leftarrow \mathcal{B}' \cup \{p\}$;
10 **return** $\mathcal{B}'$;

## C Backtracking and stitching process

Algorithm 5 shows the backtracking and stitching process to collect forward completion cost of paths. It takes an instance $(\mathcal{G}, \mathcal{R})$, the paths at each level of the beam $\{\mathcal{B}_j\}_j$ and produces, for each path in $\cup_{j=1}^J \mathcal{B}_j$, a forward completion cost.
It begins with the last beam $\mathcal{B}_j$ and iterates backwards, down to $\mathcal{B}_1$. The forward completion cost of a path $p$ at level $j$ can be obtained in two ways: by appending the completion of the path given by Rollout($p$) (line Line 4) and by appending the completion of a path at level $j + 1$ through the stitching (line Line 5-Line 9). If more than one path can be appended to form a feasible $s-t$ path, the path with the lowest cost is chosen.

<caption>**Algorithm 5: Backtrack((G, R), {$\mathcal{B}_j$})**</caption>

**Input**: $(\mathcal{G}, \mathcal{R}), \{\mathcal{B}_j\}$
**Output**: the paths in $\cup_{j=1}^J \mathcal{B}_j$ with their completion cost
1 &nbsp; $\mathcal{H} \leftarrow \emptyset, \mathcal{F}_j \leftarrow \emptyset, \forall j \in \{1, \dots, J\}$;
2 &nbsp; **for** $j \in \{J, J-1, \dots, 1\}$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp; **for** $p \in \mathcal{B}_j$ **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $p_{\text{compl}} \leftarrow \text{Rollout}(p)$;
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $a \in \mathcal{A}^+(v^+(p))$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $p' \in \mathcal{F}_{j+1} \cup \{p_t\} \wedge v^-(p') = v^+(a)$ **do**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{lsFeasible}(\mathcal{G}, \mathcal{R}, p \oplus (a \oplus p'))$ **then**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $c((a \oplus p') < c(p_{\text{compl}})$ **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $p_{\text{compl}} \leftarrow (a \oplus p'$;
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $p_{\text{compl}} \neq -$ **then**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{F}_j \leftarrow \mathcal{F}_j \cup \{p_{\text{compl}}\}$;
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{H} \leftarrow \mathcal{H} \cup \{(p, c(p_{\text{compl}}))\}$;
13 &nbsp; **return** $\mathcal{H}$

$\triangleright$ return - if infeasible path is found by Rollout
$\triangleright$ Stitching
$\triangleright$ $p_t$: empty path originated from $t$

# D Path-based feature computation

This section shows the computation of three quantities that are used for the features of the instance-path pair $\ell$ of a path $p$.

The quantity $\varphi^k(p)$ is a lower bound on the cost that completes $p$ to an $s-t$ path satisfying constraint $k$. This quantity, divided by the layer, is computed for feature (PF6). A tight evaluation of this quantity would be hard to compute, since it involves solving many resource-constrained shortest path problems. To compute a looser evaluation, we compute the minimal cost $\alpha_+^k(v)$ (resp. $\alpha_-^k(v)$) of a sequence of consecutive arcs from $v$ to $t$ in $\mathcal{G}$ having a strictly positive (resp. negative) accumulated consumption of resource $k \in \mathcal{R}$. This gives a lower bound of the cost to repair the current path if its consumption of resource $k$ does not satisfy the constraint $k$. We use the following evaluation:

$$ \varphi^k(p) = \begin{cases} \underline{h}^0(v^+(p)) & \text{if } (k \in \mathcal{R}^- \wedge r^k(p) = R^k) \vee (k \in \mathcal{R}^\le \wedge r^k(p) \le R^k), \\ \alpha_-^k(v^+(p)) & \text{if } (k \in \mathcal{R}^- \wedge r^k(p) > R^k) \vee (k \in \mathcal{R}^\le \wedge r^k(p) > R^k), \\ \alpha_+^k(v^+(p)) & \text{if } k \in \mathcal{R}^- \wedge r^k(p) < R^k, \end{cases} $$

with $\alpha_+^k(t)$ and $\alpha_-^k(t)$ computed recursively as follows:

$$ \alpha_+^k(t) = +\infty \qquad \forall k \in \mathcal{R}, $$

$$ \alpha_+^k(u) = \min \left( \min_{a \in \mathcal{A}^+(u): r_a^k > 0} \{c_a + \underline{h}^0(v^+(a))\}, \min_{a \in \mathcal{A}^+(u): r_a^k = 0} \{c_a + \alpha_+^k(v^+(a))\} \right), \forall k \in \mathcal{R}, \forall u \neq t $$

$$ \alpha_-^k(t) = +\infty \qquad \forall k \in \mathcal{R}, $$
$$ \alpha_{-}^{k}(u) = \min \left( \min_{a \in \mathcal{A}^{+}(u): r_{a}^{k} < 0} \{c_{a} + \underline{h}^{0}(v^{+}(a))\}, \quad \min_{a \in \mathcal{A}^{+}(u): r_{a}^{k} = 0} \{c_{a} + \alpha_{-}^{k}(v^{+}(a))\} \right), \forall k \in \mathcal{R}, \forall u \neq t. $$

The quantities $\kappa_{\min}^{k}(p)$ and $\kappa_{\max}^{k}(p)$ are the minimal and maximal arc cost, respectively, when attempting to satisfy constraint $k$. They are used as features (PF7) and (PF8). To calculate them, we first define the following quantities. Let $J$ be the number of topological layers of $\mathcal{G}$, $\mathcal{V}_j$ be the set of vertices of the layer $j \in \{1, \dots, J\}$ and $j(v)$ be the mapping from a vertex $v$ to the corresponding layer $j(v)$. We denote, for each topological layer $j \in \{1, \dots, J\}$, four quantities: $\mu_{\min}^{k,+}(j)$, $\mu_{\max}^{k,+}(j)$, $\mu_{\min}^{k,-}(j)$, and $\mu_{\max}^{k,-}(j)$. They represent the minimum and maximum cost of arcs that contribute positively or negatively to resource $k$, respectively, after layer $j$. These four quantities are calculated recursively as follows:

$$ \begin{aligned} \mu_{\min}^{k,+}(J) &= +\infty & \forall k \in \mathcal{R}, \\ \mu_{\min}^{k,+}(j) &= \min \left( \min_{u \in \mathcal{V}_j} \{c_a\}, \quad \mu_{\min}^{k,+}(j+1) \right), \forall k \in \mathcal{R}, \forall j \in \{J-1, \dots, 1\} \end{aligned} $$

$$ \begin{aligned} \mu_{\max}^{k,+}(J) &= -\infty & \forall k \in \mathcal{R}, \\ \mu_{\max}^{k,+}(j) &= \max \left( \max_{u \in \mathcal{V}_j} \{c_a\}, \quad \mu_{\max}^{k,+}(j+1) \right), \forall k \in \mathcal{R}, \forall j \in \{J-1, \dots, 1\} \end{aligned} $$

$$ \begin{aligned} \mu_{\min}^{k,-}(J) &= +\infty & \forall k \in \mathcal{R}, \\ \mu_{\min}^{k,+}(j) &= \min \left( \min_{u \in \mathcal{V}_j} \{c_a\}, \quad \mu_{\min}^{k,-}(j+1) \right), \forall k \in \mathcal{R}, \forall j \in \{J-1, \dots, 1\} \end{aligned} $$

$$ \begin{aligned} \mu_{\max}^{k,-}(J) &= -\infty & \forall k \in \mathcal{R}, \\ \mu_{\max}^{k,+}(j) &= \max \left( \max_{u \in \mathcal{V}_j} \{c_a\}, \quad \mu_{\max}^{k,-}(j+1) \right), \forall k \in \mathcal{R}, \forall j \in \{J-1, \dots, 1\} \end{aligned} $$

The quantities $\kappa_{\min}^{k}(p)$ and $\kappa_{\max}^{k}(p)$ are computed as follows:

$$ \kappa_{\min}^{k}(p) = \begin{cases} 0 & (k \in \mathcal{R}^{-} \wedge r^{k}(p) = R^{k}) \vee (k \in \mathcal{R}^{\le} \wedge r^{k}(p) \le R^{k}) \\ \mu_{\min}^{k,-}(j(v^{+}(p))) & (k \in \mathcal{R}^{-} \wedge r^{k}(p) > R^{k}) \vee (k \in \mathcal{R}^{\le} \wedge r^{k}(p) > R^{k}) \\ \mu_{\min}^{k,+}(j(v^{+}(p))) & k \in \mathcal{R}^{-} \wedge r^{k}(p) < R^{k} \end{cases} $$

$$ \kappa_{\max}^{k}(p) = \begin{cases} 0 & (k \in \mathcal{R}^{-} \wedge r^{k}(p) = R^{k}) \vee (k \in \mathcal{R}^{\le} \wedge r^{k}(p) \le R^{k}) \\ \mu_{\max}^{k,-}(j(v^{+}(p))) & (k \in \mathcal{R}^{-} \wedge r^{k}(p) > R^{k}) \vee (k \in \mathcal{R}^{\le} \wedge r^{k}(p) > R^{k}) \\ \mu_{\max}^{k,+}(j(v^{+}(p))) & k \in \mathcal{R}^{-} \wedge r^{k}(p) < R^{k} \end{cases} $$

## E DPRC formulation

Here we present the DPRC formulation of $1 \parallel \sum_i w_i T_i$ and TKP.
## DPRC(1||$\sum_i w_i T_i$) formulation

A dynamic program formulation of $1||\sum_i w_i T_i$ (see Section 5.1) starts to assign the last jobs to perform. Let $t_{max} = \sum_{j \in \mathcal{J}} p_j$ be the total processing time of all jobs. Let $\mathbf{y} \in \{0, 1\}^n$ be a binary vector of size $n$, $y_j$ be the $j$-th element of $\mathbf{y}$ which equals to 1 if job $j$ is performed and 0 otherwise. A state of the DP can be expressed by $(t, \mathbf{y})$, where $t = \sum_{j \in \mathcal{J}: y_j=1} p_j$ is the current time after assigning all jobs $j$ whose $y_j = 1$. The cost function $T$ can be expressed in a backward fashion:

$$ 
\begin{array}{l}
T(0, \mathbf{0}) = 0 \\
T(t, \mathbf{y}) = \min_{j \in \mathcal{J}: y_j=1} T(t - p_j, \mathbf{y} - \mathbf{1}_j) + w_j \max(0, t_{max} - t + p_j - d_j) \quad \mathbf{y} \in \{0, 1\}^n, t \in [1, t_{max}],
\end{array}
$$

where $\mathbf{1}_j$ is the indicator vector of set $\{j\}$. The minimal total tardiness is given by $T(t_{max}, \mathbf{1})$.

To formulate the problem as an RCSPP-DAG, we use state-space relaxation that considers only the state variable: time $t$, and we set $\mathcal{G}$ as the state-transition graph of this backward recurrence:

$$ 
\begin{array}{l}
T(0) = 0 \\
T(t) = \min_{j \in \mathcal{J}} T(t - p_j) + w_j \max(0, t_{max} - t + p_j - d_j) \quad t \in [1, t_{max}]
\end{array}
$$

The state space relaxation does not take into account the jobs already scheduled. A vertex $v$ can be identified by $t$, i.e., $\tau(v) = t$. Each vertex has at most $n$ outgoing arcs, each corresponding to a scheduling of a job. As a result, the number of vertices in $\mathcal{G}$ is $O(t_{max})$ and the number of arcs is $O(nt_{max})$.

The set $\mathcal{R}$ is built by a single family of $n$ equality constraints, each with 1 on the right-hand side and associated with a distinct job. To all arcs $a \in \mathcal{A}$ corresponding to the scheduling of job $j$, we associate a resource consumption $r_a^j = 1$ and $r_a^{j'} = 0, \forall j' \neq j$. The span of the constraint associated with $j$ is $[0, t_{max}]$.

Since each vertex has at most $n$ outgoing arcs, the number of candidate paths after each beam search expansion increases by a factor of at most $n$.

## DPRC(TKP) formulation

For TKP (see Section 5.2), we use the event-based dynamic program proposed by Clautiaux, Detienne, and Guillot (2021). Let $(1, \dots, 2n+1)$ be an ordered list of event indices. For each item, the list contains two events $e^{in}(i), e^{out}(i)$ that correspond to the start and end of its time window, respectively. An additional dummy event is included at position $2n+1$. Events are divided into two categories: those corresponding to the start of a time window, forming the set $\mathcal{E}^{in}$, and those corresponding to the end, forming the set $\mathcal{E}^{out}$. For each event $e \in \{1, \dots, 2n\}$, we denote by $i(e) \in \mathcal{I}$ the associated item, and by $t(e)$ the time at which the event occurs, i.e., $t(e) = s_{i(e)}$ if $e \in \mathcal{E}^{in}$, and $t(e) = f_{i(e)}$ if $e \in \mathcal{E}^{out}$. Events related to the items are ordered from 1 to $2n$ as follows: $e < e'$ if $t(e) < t(e')$ or $(t(e) = t(e')) \wedge e \in \mathcal{E}^{out} \wedge e' \in \mathcal{E}^{in}$.

To define the dynamic program, we define a state as a tuple $(e, w, \mathbf{d})$, where $e \in \{1, \dots, 2n+1\}$ is the current event, $w \in \mathbb{Z}_+$ the current knapsack capacity consumption, and $\mathbf{d} \in \{0, 1\}^n$ the characteristic vector of the set of items currently in the knapsack. We also define $\Delta_i \in \{0, 1\}^n$ as the indicator vector corresponding to item $i$, representing the inclusion of item $i$ in the knapsack. The cost function $P$ can be expressed in a backward fashion:

$$ 
\begin{array}{l}
P(2n + 1, 0, \mathbf{0}) = 0, \\
P(e, w, \mathbf{d}) = 
\begin{cases}
\min \{ P(e + 1, w, \mathbf{d}), -\frac{1}{2}p_{i(e)} + P(e + 1, w + w_{i(e)}, \mathbf{d} + \Delta_{i(e)}) \} & \text{if } e \in \mathcal{E}^{\text{in}} \wedge w + w_{i(e)} \le W, \\
P(e + 1, w, \mathbf{d}) & \text{if } e \in \mathcal{E}^{\text{in}} \wedge w + w_{i(e)} > W, \\
\min \{ P(e + 1, w, \mathbf{d}), -\frac{1}{2}p_{i(e)} + P(e + 1, w - w_{i(e)}, \mathbf{d} - \Delta_{i(e)}) \} & \text{if } e \in \mathcal{E}^{\text{out}} \wedge \mathbf{d}_{i(e)} = 1, \\
P(e + 1, w, \mathbf{d}) & \text{if } e \in \mathcal{E}^{\text{out}} \wedge \mathbf{d}_{i(e)} = 0,
\end{cases}
\end{array}
$$
The graph $\mathcal{G}$ is constructed from the state-space relaxation that only considers two state variables: the current event $e$ (the first state variable), and the current consumption $w$ of the capacity (the second variable $d=2$). The backward recurrence of the construction is as follows:

$$ P(2n+1, 0) = 0, \\ P(e, w) = \begin{cases} \min\{P(e+1, w), -\frac{1}{2}p_{i(e)} + P(e+1, w+w_{i(e)})\} & \text{if } e \in \mathcal{E}^{\text{in}} \wedge w+w_{i(e)} \le W, \\ P(e+1, w) & \text{if } e \in \mathcal{E}^{\text{in}} \wedge w+w_{i(e)} > W, \\ \min\{P(e+1, w), -\frac{1}{2}p_{i(e)} + P(e+1, w-w_{i(e)})\} & \text{if } e \in \mathcal{E}^{\text{out}} \end{cases} $$

A vertex $v \in \mathcal{V}$ is identified by an event $e$ and a consumption $w$ of capacity, i.e., $\tau(v) = e$, $\delta^2(v) = w$. Given an event $e$, all vertices associated with $e$ have at most two outgoing arcs. If $e \in \mathcal{E}^{\text{in}}$ and $w + w_{i(e)} \le W$, the two outgoing arcs correspond to the addition and the non-addition of item $i(e)$. If $e \in \mathcal{E}^{\text{in}}$ and $w + w_{i(e)} > W$, the vertex $v$ has only one outgoing arc corresponding to the non-addition of item $i(e)$. If $e \in \mathcal{E}^{\text{out}}$, the two outgoing arcs correspond to the removal and the non-removal of item $i(e)$. As a result, the number of vertices in $\mathcal{G}$ is $O(2nW)$ and the number of arcs is $O(4nW)$.

The set $\mathcal{R}$ is built by a single family of $n$ equality constraints, each with right-hand side 0 and associated with a distinct item. Arcs corresponding to the addition of item $i$ carry a consumption $r_a^i = 1$ for that resource, while arcs representing the removal of item $i$ carry a consumption $r_a^i = -1$. For resource consumption of constraint $j \neq i$ on arc $a$, we have $r_a^j = 0, \forall j \neq i$. The span of the constraint associated with item $i$ is $[e^{\text{in}}(i), e^{\text{out}}(i)]$.

Since each vertex has at most two outgoing arcs, the number of candidate paths after each beam search expansion is doubled.

## F Iterated local search for $1 \parallel \sum_i w_i T_i$

Our ILS is initialized with a solution generated by a greedy algorithm using either Earliest Due Date (EDD) rule or the Shortest Processing Time (SPT) rule. The best solution among the two runs is returned.

The ILS algorithm alternates between two phases: a local search phase that improves the current solution, and a perturbation mechanism that helps escape the local optimum. In our implementation, the local search is a variable neighborhood descent (VND). It relies on an ordered list of local search operators (neighborhoods). If an improving move is found using the current operator, the search restarts from the first operator in the list. If no improvement is found, it proceeds to the next operator. The search terminates at a local optimum when the last operator fails to improve the solution. Our VND uses four operators. The first operator performs adjacent pairwise swaps (i.e., it swaps jobs at positions $i$ and $i+1$). The second operator attempts to move a job forward by at least two positions (from position $i$ to $j \ge i+2$). The third operator tries to move a job backward by up to two positions (to $j = i-1$ or $i-2$). The fourth operator considers all possible pairwise swaps ($i < j$). When the VND cannot further improve the solution, a diversification step randomly swaps pairs of jobs, after which the VND restarts from this modified solution. The algorithm ends after reaching the maximum number of iterations, returning the cost of the best sequence found.

## G Impact of lower and upper bound scoring functions on ABSM

This section presents the results of ABSM without the integration of ML. For each of the four cases, we vary the proportion of paths (rLB) selected by the Lagrangian lower bound in $\{0, 0.1, \dots, 0.9, 1.0\}$, while setting the proportion of paths (rRoll) selected by the rollout-based upper
bound as 1-rLB. For each instance size, we report the mean and standard deviation of the gap in Tables 12, 13, 10, and 11.

<table><caption>Table 10: DPRC(1|| $\sum_i w_i T_i$) with $\lambda = \lambda^*$: tuning ratios rLB/rRoll for ABSM, best combination 0.8/0.2</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 100<br>Gap (%)</th><th colspan="2">n = 200<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>1.00</td><td>-</td><td>1.211</td><td>2.822</td><td>1.585</td><td>2.615</td></tr><tr><td>0.10</td><td>0.90</td><td>-</td><td>0.408</td><td>1.521</td><td>1.026</td><td>2.040</td></tr><tr><td>0.20</td><td>0.80</td><td>-</td><td>0.281</td><td>1.440</td><td>0.923</td><td>1.874</td></tr><tr><td>0.30</td><td>0.70</td><td>-</td><td>0.272</td><td>1.455</td><td>0.919</td><td>1.896</td></tr><tr><td>0.40</td><td>0.60</td><td>-</td><td><strong>0.250</strong></td><td>1.444</td><td>0.860</td><td>1.862</td></tr><tr><td>0.50</td><td>0.50</td><td>-</td><td>0.252</td><td>1.443</td><td>0.878</td><td>1.879</td></tr><tr><td>0.60</td><td>0.40</td><td>-</td><td>0.256</td><td>1.449</td><td>0.870</td><td>1.875</td></tr><tr><td>0.70</td><td>0.30</td><td>-</td><td>0.257</td><td>1.449</td><td>0.770</td><td>1.598</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td>0.264</td><td>1.465</td><td><strong>0.747</strong></td><td>1.576</td></tr><tr><td>0.90</td><td>0.10</td><td>-</td><td>0.300</td><td>1.576</td><td>0.799</td><td>1.642</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>0.420</td><td>1.837</td><td>0.952</td><td>2.028</td></tr></tbody></table>

<table><caption>Table 11: DPRC(1|| $\sum_i w_i T_i$) with $\lambda = 0$: tuning ratios rLB/rRoll for ABSM, best combination 0.1/0.9</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 100<br>Gap (%)</th><th colspan="2">n = 200<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>1.00</td><td>-</td><td>5.334</td><td>5.875</td><td>8.558</td><td>7.283</td></tr><tr><td>0.10</td><td>0.90</td><td>-</td><td>5.003</td><td>5.003</td><td><strong>7.734</strong></td><td>6.481</td></tr><tr><td>0.20</td><td>0.80</td><td>-</td><td>5.034</td><td>5.015</td><td>7.804</td><td>6.554</td></tr><tr><td>0.30</td><td>0.70</td><td>-</td><td>5.029</td><td>4.995</td><td>7.781</td><td>6.513</td></tr><tr><td>0.40</td><td>0.60</td><td>-</td><td><strong>4.983</strong></td><td>4.920</td><td>7.824</td><td>6.561</td></tr><tr><td>0.50</td><td>0.50</td><td>-</td><td>5.009</td><td>4.875</td><td>7.776</td><td>6.469</td></tr><tr><td>0.60</td><td>0.40</td><td>-</td><td>5.073</td><td>4.887</td><td>7.917</td><td>6.581</td></tr><tr><td>0.70</td><td>0.30</td><td>-</td><td>5.197</td><td>4.953</td><td>7.999</td><td>6.536</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td>5.354</td><td>4.988</td><td>8.212</td><td>6.638</td></tr><tr><td>0.90</td><td>0.10</td><td>-</td><td>5.689</td><td>5.237</td><td>8.566</td><td>6.743</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>13.125</td><td>11.274</td><td>26.915</td><td>20.622</td></tr></tbody></table>
<table><caption>Table 12: DPRC(TKP) with $\lambda = \lambda^*$: tuning ratios rLB/rRoll for ABSM, best combination 1.0/0.0</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 1000<br>Gap (%)</th><th colspan="2">n = 2000<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>1.00</td><td>-</td><td>0.172</td><td>0.191</td><td>0.345</td><td>0.173</td></tr><tr><td>0.10</td><td>0.90</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>0.035</td><td>0.046</td></tr><tr><td>0.20</td><td>0.80</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>0.012</td><td>0.036</td></tr><tr><td>0.30</td><td>0.70</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.002</td><td>0.027</td></tr><tr><td>0.40</td><td>0.60</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.005</td><td>0.026</td></tr><tr><td>0.50</td><td>0.50</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.007</td><td>0.025</td></tr><tr><td>0.60</td><td>0.40</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.008</td><td>0.024</td></tr><tr><td>0.70</td><td>0.30</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.009</td><td>0.023</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.009</td><td>0.023</td></tr><tr><td>0.90</td><td>0.10</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.009</td><td>0.024</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td><strong>-0.001</strong></td><td>0.005</td><td><strong>-0.011</strong></td><td>0.022</td></tr></tbody></table>

<table><caption>Table 13: DPRC(TKP) with $\lambda = 0$: tuning ratios rLB/rRoll for ABSM, best combination 0.8/0.2</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 1000<br>Gap (%)</th><th colspan="2">n = 2000<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>1.00</td><td>-</td><td>5.442</td><td>3.765</td><td>15.181</td><td>2.813</td></tr><tr><td>0.10</td><td>0.90</td><td>-</td><td>1.257</td><td>1.045</td><td>7.252</td><td>1.532</td></tr><tr><td>0.20</td><td>0.80</td><td>-</td><td>1.033</td><td>0.900</td><td>6.492</td><td>1.431</td></tr><tr><td>0.30</td><td>0.70</td><td>-</td><td>0.938</td><td>0.883</td><td>6.085</td><td>1.452</td></tr><tr><td>0.40</td><td>0.60</td><td>-</td><td>0.898</td><td>0.857</td><td>5.860</td><td>1.367</td></tr><tr><td>0.50</td><td>0.50</td><td>-</td><td>0.852</td><td>0.861</td><td>5.720</td><td>1.408</td></tr><tr><td>0.60</td><td>0.40</td><td>-</td><td>0.816</td><td>0.837</td><td>5.701</td><td>1.403</td></tr><tr><td>0.70</td><td>0.30</td><td>-</td><td>0.802</td><td>0.838</td><td>5.523</td><td>1.445</td></tr><tr><td>0.80</td><td>0.20</td><td>-</td><td>0.796</td><td>0.871</td><td><strong>5.512</strong></td><td>1.440</td></tr><tr><td>0.90</td><td>0.10</td><td>-</td><td><strong>0.789</strong></td><td>0.844</td><td>5.555</td><td>1.581</td></tr><tr><td>1.00</td><td>-</td><td>-</td><td>0.847</td><td>0.899</td><td>5.845</td><td>1.591</td></tr></tbody></table>

# H Tuning the ML model and the proportion of paths selected by ML

## H.1 Preliminary model settings

Table 14 and 15 show the preliminary feature set and model size settings for DPRC($1|| \sum_i w_i T_i$) and DPRC(TKP), each with two Lagrangian settings: whether the optimal Lagrangian multipliers are available or not. These settings are used in Section 6.2.

## H.2 Results using the preliminary models

In Tables 16, 17, 18, and 19, we report the results for four configurations of LB/Roll/ML, each presented in a separate row. The four rows show the results when ML is incorporated. We vary the proportion of the beam assigned to ML as $[0.25, 0.5, 0.75, 1.0] \times \beta$. The remaining beam size is shared between LB and Roll with the best-performing combination.
<table><caption>Table 14: Preliminary features used for DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>) and DPRC(TKP) under the two Lagrangian multiplier settings λ = λ* and λ = 0. The term "shared" in the λ column indicates that the features are used in both Lagrangian multiplier settings for a given DPRC(P). The last column stands for the number of clusters.</caption><thead><tr><th>DPRC(P)</th><th>λ</th><th>Type</th><th>Aggregators</th><th>clustering</th><th>q</th></tr></thead><tbody><tr><td rowspan="10">1||<span style="text-decoration: overline;">∑</span><sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td rowspan="3">λ*, 0</td><td>GF5 GF6 VF1 PF1 PF3</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PF8</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>GF3</td><td>skew kurt</td><td></td><td></td></tr><tr><td rowspan="6">λ*</td><td>GF3</td><td>mean std</td><td>all</td><td>-</td></tr><tr><td>GF1 GF2 GF4</td><td>mean std skew kurt</td><td>all</td><td>-</td></tr><tr><td>GF4</td><td>mean std skew kurt</td><td>q-quant-lag</td><td>4</td></tr><tr><td>PF5</td><td>mean</td><td>q-quant-lag</td><td>4</td></tr><tr><td>GF7 GF8 PF2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PF5 PF7</td><td>mean</td><td>all</td><td>-</td></tr><tr><td rowspan="10">TKP</td><td rowspan="3">λ*, 0</td><td>PF6</td><td>max</td><td>all</td><td>-</td></tr><tr><td>PF7</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>PF6</td><td>mean</td><td>all</td><td>-</td></tr><tr><td rowspan="2">λ*</td><td>GF1</td><td>mean std skew kurt</td><td>all</td><td>-</td></tr><tr><td>PF4</td><td>mean-in-span</td><td>q-quant-lag</td><td>4</td></tr><tr><td rowspan="3">λ*</td><td>GF5</td><td>-</td><td>-</td><td>-</td></tr><tr><td>GF2</td><td>mean kurt</td><td>all</td><td>-</td></tr><tr><td>VF1 VF2(2) PF1 PF3</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">TKP</td><td rowspan="2">λ*</td><td>PF4</td><td>mean-in-span</td><td>all</td><td>-</td></tr><tr><td>PF4</td><td>mean-in-span</td><td>q-quant-spanL</td><td>4</td></tr><tr><td>PF6</td><td>max</td><td>q-quant-mincost</td><td></td></tr><tr><td>PF7</td><td>mean</td><td>q-max-disjoint</td><td></td></tr><tr><td>0</td><td>PF5</td><td>mean</td><td>q-quant-maxcost</td><td>4</td></tr><tr><td rowspan="10">TKP</td><td rowspan="3">λ*, 0</td><td>GF5</td><td>-</td><td>-</td><td>-</td></tr><tr><td>GF2</td><td>mean kurt</td><td>all</td><td>-</td></tr><tr><td>VF1 VF2(2) PF1 PF3</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="2">λ*</td><td>PF4</td><td>mean-in-span</td><td>all</td><td>-</td></tr><tr><td>PF4</td><td>mean-in-span</td><td>q-quant-spanL</td><td>4</td></tr><tr><td>PF6</td><td>max</td><td>q-quant-mincost</td><td></td></tr><tr><td>PF7</td><td>mean</td><td>q-max-disjoint</td><td></td></tr><tr><td>PF6</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>PF7</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>λ*</td><td>GF1</td><td>mean std skew kurt</td><td>all</td><td>-</td></tr><tr><td>λ*</td><td>PF4</td><td>mean-in-span</td><td>q-quant-lag</td><td>4</td></tr><tr><td>0</td><td>GF5</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>Table 15: Preliminary model size</caption><thead><tr><th>DPRC(P)</th><th>λ</th><th># hid.</th><th># neu.</th></tr></thead><tbody><tr><td>1||<span style="text-decoration: overline;">∑</span><sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td>λ*</td><td>4</td><td>100</td></tr><tr><td>1||<span style="text-decoration: overline;">∑</span><sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td>0</td><td>3</td><td>50</td></tr><tr><td>TKP</td><td>λ*</td><td>3</td><td>200</td></tr><tr><td>TKP</td><td>0</td><td>3</td><td>200</td></tr></tbody></table>

<table><caption>Table 16: DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>) with λ = λ*: tuning ratio rML while the ratio rLB/rRoll is fixed according to the best combination 0.8/0.2 previously determined. Best configuration rLB/rRoll/rML = 0.2/0.05/0.75.</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 100</th><th colspan="2">n = 200</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>1.00</td><td>0.289</td><td>1.687</td><td>0.923</td><td>1.797</td></tr><tr><td>0.20</td><td>0.05</td><td>0.75</td><td><strong>0.167</strong></td><td>1.389</td><td><strong>0.474</strong></td><td>1.330</td></tr><tr><td>0.40</td><td>0.10</td><td>0.50</td><td>0.185</td><td>1.425</td><td>0.484</td><td>1.339</td></tr><tr><td>0.60</td><td>0.15</td><td>0.25</td><td>0.186</td><td>1.418</td><td>0.528</td><td>1.383</td></tr></tbody></table>

<table><caption>Table 17: DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>) with λ = 0: tuning ratio rML while the ratio rLB/rRoll is fixed according to the best combination 0.1/0.9 previously determined. Best configuration rLB/rRoll/rML = 0.05/0.45/0.5.</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 100</th><th colspan="2">n = 200</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>1.000</td><td>4.951</td><td>6.808</td><td>14.017</td><td>18.043</td></tr><tr><td>0.025</td><td>0.225</td><td>0.750</td><td>3.637</td><td>5.063</td><td>6.264</td><td>7.277</td></tr><tr><td>0.050</td><td>0.450</td><td>0.500</td><td><strong>3.586</strong></td><td>4.972</td><td><strong>6.165</strong></td><td>7.151</td></tr><tr><td>0.075</td><td>0.675</td><td>0.250</td><td>3.593</td><td>4.825</td><td>6.194</td><td>7.022</td></tr></tbody></table>
<table><caption>Table 18: DPRC(TKP) with $\lambda = \lambda^*$: tuning ratio rML while the ratio rLB/rRoll is fixed according to the best combination 1.0/0.0 previously determined. Best configuration rLB/rRoll/rML = 0.75/0/0.25.</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 1000<br>Gap (%)</th><th colspan="2">n = 2000<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>1.00</td><td><strong>-0.001</strong></td><td>0.006</td><td>0.097</td><td>0.113</td></tr><tr><td>0.25</td><td>-</td><td>0.75</td><td><strong>-0.001</strong></td><td>0.005</td><td>0.005</td><td>0.027</td></tr><tr><td>0.50</td><td>-</td><td>0.50</td><td><strong>-0.001</strong></td><td>0.005</td><td>-0.002</td><td>0.026</td></tr><tr><td>0.75</td><td>-</td><td>0.25</td><td><strong>-0.001</strong></td><td>0.005</td><td><strong>-0.008</strong></td><td>0.023</td></tr></tbody></table>

<table><caption>Table 19: DPRC(TKP) with $\lambda = 0$: tuning ratio rML while the ratio rLB/rRoll is fixed according to the best combination 0.8/0.2 previously determined. Best configuration rLB/rRoll/rML = 0.6/0.15/0.25.</caption><thead><tr><th rowspan="2">rLB</th><th rowspan="2">rRoll</th><th rowspan="2">rML</th><th colspan="2">n = 1000<br>Gap (%)</th><th colspan="2">n = 2000<br>Gap (%)</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>1.00</td><td>2.060</td><td>1.598</td><td>10.498</td><td>1.760</td></tr><tr><td>0.20</td><td>0.05</td><td>0.75</td><td>1.167</td><td>1.004</td><td>7.075</td><td>1.383</td></tr><tr><td>0.40</td><td>0.10</td><td>0.50</td><td>0.997</td><td>0.887</td><td>6.050</td><td>1.454</td></tr><tr><td>0.60</td><td>0.15</td><td>0.25</td><td><strong>0.917</strong></td><td>0.840</td><td><strong>5.824</strong></td><td>1.520</td></tr></tbody></table>
## H.3 Feature selection

This section provides the detailed results of the experiments described in Section 6.3. Table 20 summarizes how the performance of ML-ABSM changes as the less relevant features, based on SHAP values, are incrementally removed. The column "Feat. cfg" denotes arbitrary identifiers for the feature sets used in training, where larger identifiers correspond to sets with more features removed. The identifier highlighted in bold corresponds to the final feature set. In this table, we also give the specific features eliminated at each step.

<table><caption>Table 20: ML-ABSM performance of ML models for feature selection</caption><thead><tr><th rowspan="2">DPRC(P)</th><th rowspan="2">λ</th><th rowspan="2">Feat. cfg</th><th rowspan="2">#Feat.</th><th rowspan="2">Feat. rm.</th><th rowspan="2">Aggregators</th><th rowspan="2">Clustering</th><th rowspan="2">q</th><th colspan="2">Gap(%)</th></tr><tr><th>mean</th><th>std</th></tr></thead><tbody><tr><td rowspan="12">TKP</td><td rowspan="12">λ*</td><td>prelim</td><td>29</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-0.012</td><td>0.030</td></tr><tr><td>*tkp1</td><td>27</td><td>GF3<br>GF2</td><td>kurt</td><td>-</td><td>-</td><td>-0.014</td><td>0.028</td></tr><tr><td>*tkp2</td><td>25</td><td>GF2<br>GF1</td><td>mean<br>kurt</td><td>-</td><td>-</td><td>-0.010</td><td>0.027</td></tr><tr><td>*tkp3</td><td>21</td><td>PF4</td><td>mean-inspan</td><td>q-max-disjoint</td><td>4</td><td>-0.013</td><td>0.030</td></tr><tr><td>*tkp4</td><td>17</td><td>PF4</td><td>mean-inspan</td><td>q-quant-lag</td><td>4</td><td>-0.012</td><td>0.025</td></tr><tr><td>*tkp5</td><td>13</td><td>PF4</td><td>mean-inspan</td><td>q-quant-mincost</td><td>4</td><td>-0.011</td><td>0.026</td></tr><tr><td>*tkp6</td><td>9</td><td>PF4</td><td>mean-inspan</td><td>q-quant-spanL</td><td>4</td><td>-0.010</td><td>0.027</td></tr><tr><td>*tkp7</td><td>7</td><td>GF1</td><td>std<br>PF6</td><td>all</td><td>-</td><td>-0.009</td><td>0.028</td></tr><tr><td>*tkp8</td><td>5</td><td>PF4<br>PF7</td><td>max<br>mean<br>mean, skew</td><td>all<br>all<br>all</td><td>-</td><td>-0.008</td><td>0.029</td></tr><tr><td>*tkp9</td><td>3</td><td>GF1</td><td>mean</td><td>-</td><td>-</td><td>-0.010</td><td>0.027</td></tr><tr><td>*tkp10</td><td>2</td><td>VF2(2)</td><td>-</td><td>-</td><td>-</td><td>-0.016</td><td>0.026</td></tr><tr><td>*tkp11</td><td>1</td><td>VF1</td><td>-</td><td>-</td><td>-</td><td>-0.016</td><td>0.026</td></tr><tr><td rowspan="3">TKP</td><td rowspan="3">0</td><td>prelim</td><td>22</td><td>-</td><td>-</td><td>-</td><td>-</td><td>5.765</td><td>1.570</td></tr><tr><td>0tkp1</td><td>20</td><td>GF2<br>PF7</td><td>kurt<br>mean</td><td>all<br>all</td><td>-</td><td>6.538</td><td>1.474</td></tr><tr><td>prelim</td><td>48</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.911</td><td>2.375</td></tr><tr><td rowspan="4">1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td rowspan="4">λ*</td><td>*wt1</td><td>45</td><td>GF8<br>GF7<br>GF2</td><td>-</td><td>-</td><td>-</td><td>1.235</td><td>2.683</td></tr><tr><td>prelim</td><td>48</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.911</td><td>2.375</td></tr><tr><td>0wt1</td><td>11</td><td>GF2</td><td>kurt</td><td>all</td><td>-</td><td>5.914</td><td>5.900</td></tr><tr><td>0wt2</td><td>9</td><td>GF2<br>GF5</td><td>skew<br>-</td><td>all<br>-</td><td>-</td><td>6.026</td><td>6.136</td></tr></tbody></table>

The final feature set is shown in 21. In total, for $1||\sum_i w_i T_i$ with $\lambda = \lambda^*$, the final feature set contains 48 features, identical to the preliminary set. For $1||\sum_i w_i T_i$ with $\lambda = 0$, the feature set is reduced from 12 to 11 features. For DPRC(TKP) with $\lambda = \lambda^*$, the feature set is reduced to a single feature, namely the Lagrangian lower bound. Finally, for DPRC(TKP) with $\lambda = 0$, the feature set remains unchanged with 22 features, the same as the preliminary set.
<table><caption>Table 21: Final features used for DPRC(1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub>) and DPRC(TKP) under the two Lagrangian multiplier settings λ* and 0. The last column stands for the number of clusters.</caption><thead><tr><th>DPRC(P)</th><th>λ</th><th>Type</th><th>Aggregators</th><th>Clustering</th><th>q</th></tr></thead><tbody><tr><td rowspan="10">1||∑<sub>i</sub> w<sub>i</sub>T<sub>i</sub></td><td rowspan="3">λ*, 0</td><td>GF5 GF6 VF1 PF1 PF3</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PF8</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>GF3</td><td>skew</td><td></td><td></td></tr><tr><td rowspan="7">λ*</td><td>GF3</td><td>mean std kurt</td><td>all</td><td>-</td></tr><tr><td>GF1 GF2 GF4</td><td>mean std skew kurt</td><td>all</td><td>-</td></tr><tr><td>GF4</td><td>mean std skew kurt</td><td>q-quant-lag</td><td>4</td></tr><tr><td>PF5</td><td>mean</td><td>q-quant-lag</td><td>4</td></tr><tr><td>GF7 GF8 PF2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PF5 PF7</td><td>mean</td><td>all</td><td>-</td></tr><tr><td>PF6</td><td>max</td><td>all</td><td>-</td></tr><tr><td rowspan="10">TKP</td><td>0</td><td>PF5</td><td>mean</td><td>q-quant-maxcost</td><td>4</td></tr><tr><td rowspan="7">λ*, 0</td><td>PF3</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="5">0</td><td>GF2</td><td>mean kurt</td><td>all</td><td>-</td></tr><tr><td>GF5 VF1 VF2(2) PF1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PF4</td><td>mean-in-span</td><td>all</td><td>-</td></tr><tr><td>PF4</td><td>mean-in-span</td><td>q-quant-spanL</td><td>4</td></tr><tr><td></td><td></td><td>q-quant-mincost</td><td></td></tr><tr><td></td><td></td><td>q-max-disjoint</td><td></td></tr><tr><td></td><td>PF6</td><td>max</td><td>all</td><td>-</td></tr><tr><td></td><td>PF7</td><td>mean</td><td>all</td><td>-</td></tr></tbody></table>