Learning to Handle Complex Constraints  for Vehicle Routing Problems

Authors: jieyi  bi
yining  ma
jianan  zhou
wen  song
zhiguang  cao
yaoxin  wu
jie  zhang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:17:34

Abstract: Vehicle Routing Problems (VRPs) can model many real-world scenarios and often involve complex constraints. While recent neural methods excel in constructing solutions based on feasibility masking, they struggle with handling complex constraints, especially when obtaining the masking itself is NP-hard. In this paper, we propose a novel Proactive Infeasibility Prevention (PIP) framework to advance the capabilities of neural methods towards more complex VRPs. Our PIP integrates the Lagrangian multiplier as a basis to enhance constraint awareness and introduces preventative infeasibility masking to proactively steer the solution construction process. Moreover, we present PIP-D, which employs an auxiliary decoder and two adaptive strategies to learn and predict these tailored masks, potentially enhancing performance while significantly reducing computational costs during training. To verify our PIP designs, we conduct extensive experiments on the highly challenging Traveling Salesman Problem with Time Window (TSPTW), and TSP with Draft Limit (TSPDL) variants under different constraint hardness levels. Notably, our PIP is generic to boost many neural methods, and exhibits both a significant reduction in infeasible rate and a substantial improvement in solution quality.

Main paper: 
 introduction: 
 Vehicle routing problems (VRPs) are NP-hard combinatorial optimization problems with complex constraints that model real-world scenarios, such as logistics [1] and supply chains [2] . For decades, traditional solvers relied on hand-crafted rules for VRP optimization and constraint handling. Recently, the learning-to-optimize community [3] has successfully trained deep neural networks to automatically construct VRP solutions in an end-to-end manner [4] [5] [6] . These data-driven neural methods offer greater efficiency and high parallelism for batch optimization, making them favorable alternatives.In general, neural methods construct VRP solutions by autoregressively sampling a node from its predicted distribution while masking out nodes that would violate constraints to ensure the solution's feasibility. Despite successes (e.g., on TSP and CVRP), this masking mechanism assumes that 1) the feasibility of the entire solution can be properly decomposed into the feasibility of each node selection † Yining Ma is the corresponding author. step, and 2) ground truth masks are easily obtainable for each step. However, such assumptions may fail in VRPs (e.g., TSPTW) with complex interdependent constraints among decision variables (i.e., nodes). As will be discussed in Section 4.1, this creates a masking dilemma -considering only the local feasibility of node selections does not guarantee the overall feasibility of the constructed solutions, while computing global feasibility masks that account for future impacts transforms masking itself into another intractable NP-hard problem.These observations highlight significant gaps in applying recent neural methods to practical VRPs, necessitating research on new constraint-handling frameworks. In the literature, few studies have focused on novel ways of handling feasibility in neural constructive solvers. Although preliminary methods have attempted to mitigate it by relaxing constraints into soft ones [7, 8] or supplementing networks with more feasibility-related features [9] , the former is prone to failure when applied to more complex scenarios, while the latter requires problem-specific features and a large supervised learning dataset, limiting its adaptability to broader VRPs. Consequently, neural methods still show limited flexibility, poor feasibility rates and large optimality gaps in solving those complex VRPs.In this paper, we propose a novel Proactive Infeasibility Prevention (PIP) framework to extend the capabilities of neural constructive methods for VRPs with complex interdependent constraints. Our PIP first integrates the Lagrangian multiplier method into the reinforcement learning framework of neural methods, promoting initial constraint awareness and search guidance. To further address the limitations of the Lagrangian multiplier method on complex constraints, we then introduce preventative infeasibility masking to proactively steer the search to (near-)feasible regions during solution construction. By doing so, PIP significantly enhances feasibility rates and reduces optimality gaps. Moreover, to reduce the costs of obtaining preventative infeasibility information during training, we present PIP-D, which employs an auxiliary decoder to learn and predict masking information. Our PIP-D also incorporates two adaptive strategies: one to balance infeasible and feasible masking information for different problem hardness, and another to periodically update the model so as to balance training efficiency with prediction accuracy. These advancements enable PIP-D to achieve comparable or even better performance than PIP, particularly on larger and more constrained VRP instances, while significantly reducing computational complexity.Our contributions are as follows: 1) Conceptually, we represent an early work to address and advance the handling of complex interdependent constraints in VRPs, where the original masking loses effectiveness due to the aforementioned dilemma, thereby extending the applications of neural methods to more practical scenarios. 2) Methodologically, we propose novel PIP and PIP-D approaches that can boost the capabilities of most constructive neural methods. Specifically, we leverage the Lagrangian multiplier method and introduce preventative infeasibility masking, which is further learned by an auxiliary decoder network with two adaptive strategies, to proactively and efficiently steer the search during solution construction. 3) Experimentally, we conduct extensive validation to demonstrate the effectiveness and versatility of PIP across various backbone models (i.e., AM [4] , POMO [5] , and GFACS [10] ) and complex VRP variants (i.e., TSPTW and TSPDL). Notably, PIP achieves both a significant (up to 93.52%) reduction in infeasible rate and a substantial improvement in solution quality on synthetic and benchmark datasets with different constraint hardness levels. 
related work: 
 Neural solvers for VRPs. Existing literature on learning to optimize VRPs features two primary paradigms: constructive solvers and iterative solvers. Constructive solvers learn policies to construct solutions from scratch in an end-to-end manner. Early works introduce Pointer Network to approximate the optimal solution to TSP [11, 12] and CVRP [13] in an autoregressive (AR) way. Among all AR solvers, the attention-based model (AM) [4] represents a milestone in solving a series of VRPs. Later, the policy optimization with multiple optima (POMO) [5] further improves upon AM by considering the symmetry property of VRP solutions. Numerous recent studies have then aimed to further enhance their performance [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] and versatility [24] [25] [26] [27] . Besides the AR methods, several works construct a heatmap, which indicates the probability distribution of each edge being part of the optimal solution, to solve VRPs in a non-autoregressive (NAR) manner [10, [28] [29] [30] [31] [32] [33] [34] . Despite the superior performance on large-scale instances, we note that a recent work [35] questions the effectiveness of heatmap generative methods due to the misalignment of training and testing objectives. Differently, iterative solvers learn policies to iteratively refine an initial solution. The policies are often trained in contexts of classic heuristics or meta-heuristics for obtaining more efficient and effective search components [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] . Generally, constructive solvers can efficiently achieve desirable performance levels, whereas iterative solvers hold the potential to search for nearoptimal solutions with a prolonged time budget. Additionally, there are also several works studying the scalability [46] [47] [48] [49] [50] [51] [52] , generalization [53] [54] [55] [56] [57] [58] , and robustness [59, 60] of neural VRP solvers, and leveraging large language models (LLMs) to optimize VRPs [61] [62] [63] .Constraint handling for VRPs. Most neural methods for VRPs manage constraints using a feasibility masking mechanism that eliminates actions leading to infeasible solutions during construction or iteration search [4, 28, 32, 41] . However, such a mechanism assumes the availability of accurate masks and often lacks constraint awareness learning during training, which is not always practical or desirable. For example, Zhao et al. [64] highlighted the benefits of learning to modulate agent behaviours in the 3D Bin Packing Problem, and Ma et al. [45] showed that temporary constraint violations could enhance neural iterative solvers. Despite their successes, these approaches are inherently unsuitable for assisting constructive solvers to address the VRPs with complex interdependent constraints studied in this paper. While Tang et al. [8] and Zhang et al. [7] proposed methods to transform hard constraints into soft ones via relaxation techniques and problem redefinition, respectively, they may only be able to yield near-feasible solutions with large infeasible rates for VRPs with complex constraints. More recently, Chen et al. [9] developed a multi-step look-ahead (MUSLA) method specifically tailored for TSPTW, incorporating problem-specific features and a large supervised learning dataset. In contrast, this paper proposes a more flexible and generic PIP framework based on novel ideas of preventative infeasibility masking, learnable decoders, and adaptive strategies to advance a broader range of neural methods without needing labelled training data. 
preliminaries: 
 In this paper, we mainly consider two VRP variants with complex interdependent constraints (i.e., TSPTW and TSPDL), and neural solvers (i.e., AM [4] , POMO [5] and GFACS [10] ).Problem definitions and notations. A VRP instance can be defined over a complete graph G = {V, E}, where V = {v 0 , v 1 , . . . , v n } denotes the node set, and E = {e (v i , v j ) |v i , v j ∈ V, i ̸ = j} denotes the directed edge set among all nodes. The objective is to minimize the total cost (e.g. Euclidean length) of the solution tour. To form a feasible solution, each node in V should be visited exactly once while respecting problem-specific constraints. We consider two types of VRP constraints that are practical in industry: 1) Time window constraint: The arrival time at node v i , denoted as t i , must fall within a customer-specific time window [l i , u i ]. If The vehicle arrives early (i.e., t i < l i ), it must wait until l i ; 2) Draft limit constraint: Each node v i represents a port with a non-negative demand δ i and a maximum draft d i . We denote the current cumulative load of the freighter at port v i as α i in a given solution, which should not exceed the corresponding maximum draft d i of the port.Constructive solvers for VRPs. Popular neural constructive solvers [4, 5] typically parameterize the policy using an encoder-decoder model with parameter θ, trained with reinforcement learning (RL). Given a VRP instance G = {V, E}, the features of each node v i are represented as f v i = {x i , y i , c i }, where x i , y i are node coordinates, c i represents constraint-related features (e.g., c i = {l i , u i } for time windows in TSPTW and c i = {δ i , d i } for demand and draft limits in TSPDL). The encoder transforms node features into high-dimensional representation embeddings h i , which, combined with the context of the partial tour, represent the current state. The decoder takes them as inputs and outputs probabilities for candidate nodes (actions). The reward R(τ |G) is the negative tour length. The policy π θ is typically trained using REINFORCE [65] as follows:EQUATIONwhere K denotes the number of sampled solutions τ i for a given training instance G, and b(•) is a baseline function to reduce the variance. Specifically, the baseline is the reward (negative tour length) of the solution derived greedily in AM or the average reward of sampled solutions 1K K i=1 R (τ i |G) in POMO.Notably, POMO stipulates the starting node of each solution for diversification, which, however, may hinder solution feasibility in our studied complex constrained problems. Based on our preliminary experiments, POMO with and without diverse starting nodes achieve around 50.70% and 1.75% infeasible rates on the easy TSPTW-50 datasets, respectively. Therefore, we remove the starting node stipulation in POMO and instead sample n solutions to calculate the baseline. Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [NA]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
methodology: 
 We now discuss the limitations of existing masking mechanisms in solving VRPs with complex interdependent constraints, followed by a detailed introduction of our PIP and PIP-D frameworks.The core of feasibility masking in neural constructive solvers is to filter out invalid actions that violate constraints, based on the assumption that the global feasibility can be decomposed into the feasibility of each node selection step, and that ground truth masks are obtainable for each step. Without loss of generality, we illustrate the dilemma of feasibility masking using a TSPTW example. In TSPTW, nodes are masked out if they have been visited or cannot be visited before their time window closes. However, the feasibility of selecting a node at a particular step impacts the current time, thereby affecting all future selections due to the interdependence of time window constraints. Thus, considering only local feasibility does not guarantee overall feasibility and may lead to irreversible infeasibility. For instance, in a 4-node TSPTW instance with time windows {[0, 7], [1, 4] , [2, 4] , [2, 6] } as illustrated in the left panel of Figure 1 , there is a feasible solution τ= (v 0 → v 1 → v 2 → v 3 ).Yet, with the partial solution v 0 → v 1 , both v 2 and v 3 appear locally feasible. If the solver selects v 3 , the tour becomes infeasible irreversibly. A potential remedy is to compute global feasibility masks that consider all future possibilities, as illustrated in the right panel of Figure 1 . However, this makes masking itself an NP-hard problem, which creates a dilemma between ensuring solution feasibility and managing computational complexity. Note that this dilemma is less critical in CVRPTW, which involves multiple vehicles and routes, providing more flexibility. If one route becomes infeasible, another vehicle departing at time 0 can cover the missed nodes, reducing the impact of constraint interdependencies. However, this issue is severe in TSPTW and other variants like TSPDL.We first formulate the solution construction process of VRP as a Constrained MDP (CMDP) defined by the tuple (S, A, P, R, C), where S is the state space, A is the action space that travels from node v i to node v j , R : S × A × S is the reward function, C : S × A × S is the constraint violation cost (penalty) function, and P : S × A × S → [0, 1] is the transition probability function. At each time step, the neural solver outputs the probability of all candidate nodes, and selects one to construct a complete solution τ . The objective of CMDP is to learn a policy π θ : S → P(A) that maximizes the summation of the state-wise reward subject to certain constraints,EQUATIONwhere J is the expected return of the policy, Π F denotes the set of all feasible policies, κ m represents the boundary of the inequality constraints C m , and M is the number of constraints. Specifically, a feasible policy π is one whose expected value of constraint violation w.r.t C m , denoted as J Cm (π), does not exceed κ m . Note that κ m is set to 0 throughout this paper since we consider the hard constraints that do not tolerate any violation. Moreover, we set the reward function R to the negative value of the Euclidean distance between two nodes, i.e., R (eΠ Π F ! !! ! !" ! !# (b) Π Π F ! !! ! !" ! !# (a) Π Π F (c) Π Π F ! !! ! !" ! !# (d) Π Π F ! !! ! !" ! !# (e) ! !! ! !" ! !#(v i , v j )) = -||v i -v j || 2 .By applying feasibility masking, the search is confined to only feasible regions, allowing neural methods to focus solely on the objective function in Eq. ( 2) without explicitly considering constraint awareness or constraint violations. However, these methods lose effectiveness when such masks are unavailable, leading to inefficient searches in large infeasible regions. To address this, we propose PIP, combining a Lagrangian multiplier for constraint awareness and preventative infeasibility masking to confine the search space to near-feasible regions for complex constrained problems.Lagrangian-assisted constraint awareness. We design a Lagrangian multiplier based method to incorporate constraints C into the reward function R. Based on the Lagrangian Multiplier Theorem, the CMDP formulation in Eq. ( 2) is transformed into the following MDP formulation for VRPs:EQUATION)where L is the Lagrangian function, and λ m is a non-negative Lagrangian multiplier. Generally, the constraint violation term is calculated as the total violation value of all constraints. In TSPTW, J TW (τ ) = n i=0 max(t i -u i , 0), and in TSPDL, J DL (τ ) = n i=0 max(α i -d i , 0). Additionally, we introduce the number of infeasible nodes in the solution τ , termed as J IN , as an extra term in the Lagrangian function for better constraint awareness, which is empirically found to be effective to reduce the infeasibility rate. While Lagrangian relaxation has been explored in neural iterative methods for soft objectives [8] , our approach introduces a distinct constraint violation cost function tailored for neural constructive methods and considers fixing the Lagrangian multiplier λ (the dual variable) and optimizing the primal variable θ, significantly reducing computational overheads.Preventative infeasibility (PI) masking. As depicted in Figure 2 (b), the customized Lagrangian multiplier guides the neural policy towards a potentially feasible and high-quality space using Eq. (3). However, for more complex cases shown in Figure 2 (c), neural solvers may still struggle to navigate the large search space. To further improve training efficiency and solution feasibility, we introduce preventative infeasibility (PI) masking to proactively avoid selecting infeasible nodes during the solution construction process. As shown in the left panel of Figure 3 , if selecting a candidate node (i.e., orange node) results in any remaining candidates (i.e., green node) becoming potentially unvisitable in the next step due to constraint violations, it is marked as infeasible (i.e., red node) since selecting it would cause irreversible future infeasibility (see Appendix A.3 for a detailed example). Note that we employ a simple yet effective one-step PI masking in this paper to balance computational costs without iterating over all future possibilities (which is NP-hard). Together with the customized Lagrangian multiplier, our PIP proactively reduces the search space to a near-feasible domain Π F , as shown in Figures 2(d)-(e) . Notably, such PIP design is generic and can be applied to enhance most neural constructive solvers for VRPs with complex interdependent constraints. Benefiting from the constraint-aware optimization guided by our PIP, neural methods gain enhanced capabilities to address complex constraints, significantly boosting feasibility and optimality. However, acquiring the above PI information introduces extra computational costs (see Section 5). To alleviate this, we propose an auxiliary decoder network to learn and predict these masks, replacing the timeconsuming process of generating PI information with a much faster forward pass of the PIP decoder. This further accelerates the training process, resulting in an enhanced version of our PIP framework, termed PIP-D. The overall framework of PIP-D is illustrated in Figure 3 .Auxiliary PIP decoder. As presented in Figure 3 , we incorporate an auxiliary decoder to learn and predict the PI masks. Our PIP-D simultaneously involves training a routing decoder (the original one) that maximizes the expected reward of solutions in Eq. ( 3) and a PIP decoder that minimizes the prediction error on the PI masking using a weighted binary cross-entropy loss. The combined loss function is a weighted sum of these two objectives, i.e., L = αL RL + βL PIP . To ensure the generality, the PIP decoder mirrors the architecture of the backbone model and only adjusts the final output layer with Sigmoid activation. More details of our PIP decoder are provided in Appendix B.PIP-D training with adaptive strategies. Nevertheless, efficiently training the PIP decoder together with the routing decoder necessitates effective designs. We address this with two adaptive strategies. Firstly, training the PIP decoder at every gradient step would result in higher computational complexity than the original PIP, counteracting our goal of reducing training complexity. Hence, we adopt a periodic update strategy that intermittently updates the PIP decoder instead of continuously doing so. This approach is based on the observation that the PI masks recommended by the neural network tend to remain robust over short training periods. Specifically, we first train the PIP decoder with E init epochs, then periodically update E u epochs per E p epochs, and finally conduct E l -epoch updates. In this way, the computational costs are reduced and can be adaptively adjusted. Secondly, we consider balancing feasible and infeasible PI signals for instances with different inherent hardness. Given that the proportion of PI signals identified for feasible and infeasible nodes can vary significantly across different VRP variants with different inherent hardness, we employ a weighted balancing strategy to mitigate the influence of label imbalance [66] , which is formulated as follows:EQUATION)where T is the total decoding step to construct a complete solution. The weights of each category are calculated by their corresponding sample number, i.e,ω infsb = Ninfsb+Nfsb 2Ninfsb , ω fsb = Ninfsb+Nfsb 2Nfsb, where N infsb and N fsb are the number of infeasible and feasible nodes identified by our PI masking (g t ) in a specific decoding step t, respectively. Moreover, beyond the above two critical strategies, we explore additional strategies to accelerate the training of the PIP decoder, including fine-tuning and early-stopping techniques, which are discussed in Appendix D.3. # Results are adopted from [9] due to unavailable source code, with our 'Easy' settings corresponding to their 'Medium' dataset. / The corresponding results are not available due to no feasible solutions or not given by [9] . † We report the average results for instances where feasible solutions were found, which vary across different models. Despite these variations, the results for overlapping feasible instances consistently show similar patterns (see Appendix D.2). 
experiments: 
 In this paper, we propose a Proactive Infeasibility Prevention (PIP) framework and its enhanced version, PIP-D, to address the limitations of existing masking mechanisms for handling complex constraints. Notably, our PIP and PIP-D are generic and can be applied to boost various problem variants and neural methods. To evaluate the effectiveness of our method, we apply our PIP frameworks to two representative AR constructive methods, AM [4] and POMO [5] , and the latest NAR constructive GFACS [10] . For the benchmark problem, we consider two representative complex VRP variants with strong interdependent constraints that challenge existing neural methods (i.e., TSPTW and TSPDL, each at varying levels of hardness) with small problem scale n = 50, 100 for AM [4] and POMO [5] and large scale n = 500 for GFACS [10] . All the experiments are conducted on servers with NVIDIA GeForce RTX 3090 GPUs and Intel(R) Xeon(R) Gold 6326 CPU at 2.90GHz. Our implementation in PyTorch are publicly available at https://github.com/jieyibi/PIP-constraint.Implementation details. We generate instances at different hardness levels following prior works.For TSPTW [7, 9, 30, 67] , we generate three types of instances: Easy, Medium and Hard, by adjusting the width and overlap of the time window. For TSPDL [68] [69] [70] , we consider two levels of hardness: Medium and Hard. More details of such instance generation are provided in Appendix A. To ensure a comprehensive comparison, we also train and evaluate the models learned solely using our designed Lagrangian multiplier method. Meanwhile, we mark the models that use the Lagrangian multiplier with an * for clarity. For our proposed approaches, our PIP models build on the Lagrangian multiplier by further incorporating one-step preventative infeasibility masking, while the enhanced version, PIP-D, is trained with a periodically and adaptively updated PIP decoder as previously described.Hyper-parameters for training follow the original settings of the backbone models except for the ones related to the added PIP decoder. Detailed hyper-parameters and additional results are available in Appendix C and D. During inference, we adhere to the settings of the original backbone models. For the AM series models, we sample 1280 solutions per instance; for the POMO series models, we use a greedy strategy with 8× augmentation; and for the GFACS series models, we employ 100 ants to generate solutions for each instance over 10 pheromone iterations.Baselines. We compare our proposed PIP framework with two types of baselines: 1) heuristic methods, including LKH3 [71] , a strong solver designed for multiple VRP variants; OR-Tools [72] , a more flexible solver allowing different combinations of multiple diverse constraints; and Greedy Heuristics that selects locally optimal candidates at each step, where Greedy-L picks the nearest candidate and Greedy-C chooses based on complex constraints: in TSPTW, the soonest time window ends relative to the current time; and in TSPDL, the minimal draft limit; 2) Neural methods, including the original AM [4] , POMO [5] and GFACS [10] , as well as JAMPR [73] , adapted by [9] to solve TSPTW from VRPTW; and MUSLA [9] , a prior work on TSPTW trained in supervised manner, where OSLA is its one-step version and MUSLA adapt adopts an adaptive inference strategy. More details on the compared baselines are presented in Appendix C.Evaluation metrics. In this paper, we report the following metrics to evaluate the performance of our proposed PIP framework: 1) the ratio of infeasible solutions (Infeasible%), which includes the solution-level (Sol.) infeasible rate that considers all generated solutions during inference and the instance-level (Inst.) infeasible rate that considers the comprehensive results of N s solutions generated by the sampling (N s = 1, 280 in AM series models) or augmentation (N s = 8 in POMO series models). If at least one feasible solution is found among these N s solutions, the instance is considered to have feasible solutions; 2) average optimality gap (Gap) w.r.t the strong baseline LKH [71] for the best feasible solutions within N s solutions; 3) average tour length (Obj.) of the feasible best solutions within N s solutions; and 4) inference time, where we report the total time taken to solve 10,000 (n = 50 and 100) or 128 (n = 500) instances, with batch parallelism enabled on a single GPU. For baselines run in CPU, we exhibit the results in parallel on 16 CPU cores.The performance comparison on TSPTW and TSPDL at various levels of problem hardness is presented in Table 1 and Table 2 , respectively. Notably, the original backbone models AM and POMO could not solve the problem even at the easiest level. By incorporating the Lagrangian multiplier (indicated by *), the models begin to generate some feasible solutions. However, this advantage diminishes under more complex constraints. For example, the instance-level infeasibility rates for POMO and AM on Hard TSPTW-100 reach 100% in Table 1 , which is dramatically reduced to 6.28% with the addition of PIP-D, while also improving solution quality. Compared to traditional heuristics like ORTools, Greedy-L, and Greedy-C, our PIP-D consistently outperforms these methods and shows favourable results against JAMPR and MUSLA, especially in large-scale problems. Furthermore, compared to PIP, our PIP-D delivers competitive or even better objective values and optimality gaps while significantly enhancing training efficiency (e.g., 1.5 times faster for n = 50 and 5.8 times faster for n = 100, w.r.t POMO* + PIP). Notably, the superiority of PIP-D is more significant on the more constrained hardness levels and larger problem sizes. For TSPDL, we observe similar patterns, where our PIP and PIP-D models consistently outperform other baselines in terms of both infeasibility reduction and solution quality. These results validate that our PIP approach significantly reduces infeasible rates and substantially improves solution quality compared to existing neural methods. We further evaluate the capability of solving largescale problems by implementing our PIP framework on GFACS [10] . As displayed in Table 3 , equipping GFACS with our PIP significantly reduces the infeasible rate, for both the solution level and instance level and simultaneously enhances solution quality. Notably, GFACS* + PIP-D almost guarantees to obtain all feasible solutions. Different from AM and POMO, GFACS is a NAR constructive solver, which showcases the generality of our framework.Ablation on each PIP and PIP-D design. We now provide in-depth discussions on the effectiveness of the three proposed designs: the Lagrangian multiplier (*), the PI masking (PIP) and the learnable decoder (PIP-D). As shown in Tables 1 and 2 , in Easy datasets, the solution-level infeasible rate for POMO* is 2.11%, improving to 0.06% with PIP-D. This shows that, for less complex constraints, the Lagrangian multiplier alone effectively guides the policy to feasible regions, hedging the impact of PIP and PIP-D, which aligns with Figure 2 (b) and (d) where Π F is relatively large compared to Π. However, in more complex scenarios, where Π F is much smaller relative to Π (as in Figure 2 (c)), the neural policy struggles even with the Lagrangian multiplier. In such cases, our PIP and PIP-D become crucial, significantly confining the search space as depicted in Figure 2 (e). In Medium datasets, the infeasible rate drops from 18.7% in POMO* to 3.34% in POMO* + PIP-D; in Hard datasets, it drops dramatically from 100% in POMO* to 6.48% in POMO* + PIP-D. These results verify that our PIP framework achieves significant improvement, especially as problem complexity increases. Ablation on the terms in Lagrangian function. In Figure 4 , we exhibit the results with and without the J IN in Eq.( 3), which validates its efficacy of enhancing the constraint awareness.Ablation on weighted balancing strategy. Recall that the ratio of infeasible to feasible samples in PIP labels varies with the inherent constraint hardness. Our preliminary experiments suggest that such a ratio can reach up to 20:1 in the case of Hard datasets, causing significant label imbalance. This imbalance may significantly impact the performance of POMO* + PIP-D on several datasets, especially the harder ones, leading to 0% prediction accuracy on the minority class and causing a 100% infeasible rate for the backbone solver without a weighted balancing strategy. This indicates that the hardness-adaptive label balance strategy is essential. Moreover, for the accuracy of PIP-D, please refer to Appendix D.4.Ablation on periodical update strategy. In Figure 5 , we evaluate PIP-D models with fewer updates. Results show that more updates improve performance, despite a slight increase in training time.Ablation on different step numbers. Instead of iterating over all future possibilities, we use one-step PI masking to approximate NP-hard feasibility mask and reduce computational cost. To provide a Comparison with LKH3 under different inference time budget. To provide a more comprehensive comparison with LKH3, we provide additional results of LKH3 with identical time limits as the proposed approach across varying instance difficulty levels (Easy, Medium, Hard) and scales (n = 50, 100 nodes). The time limits are configured in two ways: matching the per instance inference time without parallelization and matching the total inference time with parallelization on a GPU. As shown in Table 6 , POMO*+PIP(-D) outperforms LKH3 on Hard datasets, while maintaining competitive results on Easy and Medium datasets. While comparing per-instance time might seem fair for CPU-based LKH3, ignoring parallelization could disadvantage GPU-based solvers (thus not fair for GPU-based solvers). To further explore this, we conduct another experiment but with a similar total inference time limit across both methods, which is a common practice in most existing NCO papers. Results, in Table 7 , show that our POMO*+PIP(-D) performs consistently better than LKH3 across most of the hardness. Moreover, to leverage the strengths of both approaches and further reveal the practical usage of our method, we explore a hybrid method that combines our PIP(-D) framework with LKH3. The results, as shown in Appendix D.3, reveal that LKH3's search efficiency can be significantly enhanced when initialized with solutions from our PIP(-D) framework. 
conclusions: 
 In this paper, we study an unsolved challenge in neural VRP solvers and correspondingly propose a novel Proactive Infeasibility Prevention (PIP) framework to advance their capabilities towards addressing VRPs with complex constraints. Technically, we introduce a Lagrangian multiplier method and preventative infeasibility masking to proactively guide the solution construction process. By further incorporating an auxiliary decoder, our PIP framework enhances training efficiency while exhibiting superior performance on more complex datasets. While our PIP is generic and has shown great ability to boost both AR and NAR constructive methods, one potential limitation is that it may not improve performance on all backbone solvers and all VRP variants. Future directions include: 1) exploring other strategies to reduce computational complexity, such as employing a trainable heatmap to confine the candidate space of PI masking calculation, 2) applying PIP to more neural methods at larger scales, 3) extending PIP to neural iterative solvers, 4) applying PIP to more VRP variants with complex constraints, including those hard-constrained VRPs whose feasibility masking is not NP-hard but with large optimality gaps, 5) exploring the applications of PIP in other domains, such as job shop scheduling, where operations need to be completed in a specific order and infeasibility can be proactively prevented using PIP, and 6) developing theoretical justifications for PIP. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: Licenses for existing assets are detailed in Appendix F.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: A Details of considered VRPsIn this paper, we mainly consider two VRP variants with complex interdependent constraints: Traveling Salesman Problem with Time Window (TSPTW) and TSP with Draft Limit (TSPDL). We first detail their corresponding data generation process, then demonstrate the interdependent nature of constraints inherent in these problem variants.Each TSPTW instance includes a depot node and n customer nodes with four properties: 2dimensional coordinates in Euclidean space (x i , y i ), lower bound (l i ) and upper bound (u i ) of time windows. In prior works [7, 9, 30, 67] , they generate the coordinates following a uniform distribution confined in a square box (x i , y i ) ∼ U[0, 100] while generating time window differently. Concretely, there are three ways to synthesize the time window: 1) construct a near-optimal TSP solution first and generate the time window according to the distance between the two adjacent nodes in the pre-generated near-optimal solution (e.g. in [74, 75] ); 2) construct random node permutation first and generate time window according to the distance between the two adjacent nodes in the pre-generated random solution (e.g. in [30, 67, 76] ); and 3) generate the time window under a uniform distribution without prior TSP permutation (e.g. in [7, 9] ). The first two methods generate the time window based on a pre-generated TSP solution which can guarantee the existence of feasible solutions for the generated instances. However, they diminish the impact of time window constraints due to the strong prior knowledge of TSP and the pre-generated solutions. Meanwhile, the first method necessitates obtaining a near-optimal solution for TSP initially, which incurs additional computational costs. In contrast to the first two methods, the third method appears more generic but does not guarantee feasibility. Below, we describe the detailed settings for instance generation, considering three different levels of hardness, as outlined in the main paper.Easy TSPTW. We mainly follow the settings from recent works [7, 9] and employ the third method presented above to generate time window. In specific, the lower bound of the time window l i follows a uniform distribution, i.e., l i ∼ U[0, T N ], where T N is an estimator of expected tour length in relation to the problem scale. For example, T 20 ≈ 10.9 [9] . The upper bound of the time window u i is generated based on l i , whereu i ∼ l i + T N • U[α, β],and α and β are set to 0.5 and 0.75, respectively.Medium TSPTW. It follows the same settings of the easy TSPTW, except that α and β are set to 0.1 and 0.2, respectively. To decrease α and β, we derive TSPTW instances with tighter time windows, resulting in an increased hardness level.Hard TSPTW. Different from easy and medium TSPTW, hard instances are generated following the settings of the benchmark dataset [67] . The second method is leveraged to generate time window. Concretely, we first obtain a permutation τ by randomly shuffling nodes. Then, the time window is generated based on τ following a uniform distribution,l i ∼ U [ψ i -η, ψ i ] and u i ∼ U [ψ i , ψ i + η],where ψ i denotes the cumulative distance of the partial solution until time step i. η is a factor to control the width of the time window. As preliminary experiments [9] suggest, the original η was set at 500, allowing the instance to achieve optimality using a simple greedy heuristic. Hence, in this paper, we set η to 50 to generate tighter time windows and increase hardness.Following the conventions [4, 5] , we normalize the node coordinates into [0, 1] by dividing a scale factor ρ = 100. Pertaining to the time window, we also divide l i and u i by ρ. In specific, following [30] , we first modify the upper bound of the depot u 0 to the maximum value of the upper bound of time window among all the customer nodes plus the travelling distance between them, i.e, u 0 = max(u i + ||v i -v 0 || 2 ), i ∈ [1, n].Then we use u 0 as the normalization factor, scaling all l i and u i by it to confine their values within the range of [0, 1].TSPDL is prevalent in marine transportation scenarios, which consider the vessel capacity of the freighters. Each node has its own demand δ i and draft limit d i . A TSPDL instance is derived from a TSP instance by mutating the draft limit to less than the total demand for σ% of nodes, i.e.,d i ∼ [δ i , n j=0 δ j ],while the draft limits of remaining nodes are equal to the total demand, i.e., d i = n j=0 δ j . We can adjust σ to manipulate the hardness of the TSPDL dataset. The demand is set to one for customer nodes or zero for depot node v 0 following benchmark datasets [68, 69, 71] . The mutation proportion σ% varies, being either a random value or fixed to 10, 25, 50, or 75 percent. We observe that a small σ% results in a simple TSPDL dataset that can be effectively managed with heuristics. In this paper, we focus on relatively hard problems by setting σ% to 75% for the Medium dataset and 90% for the Hard dataset. Note that the availability of feasible solutions can be guaranteed through the feasibility check, as demonstrated in the following PyTorch implementation. node_demand = torch.cat([torch.zeros((batch_size, 1)), torch.ones((batch_size, problem_size -1))], dim=1) # shape: (batch_size, problem_size) demand_sum = node_demand.sum(dim=1).unsqueeze (1) for i in range(batch_size): feasible = False while not feasible: # mutation of dl randomly occurs in w% of the nodes except the depot mutation = torch.randint (1, demand_sum[i] .int().item(), size = (problem_size * sigma // 100,)) count = torch.bincount(mutation) count_cumsum = torch.cumsum(count, dim=0) feasible = (count_cumsum <= torch.arange(0, count.size(0))).all() To summarize, the feasibility guarantee based on the data generation rules described above, and its corresponding accessibility by LKH3 [71] are presented in Table 8 . For example, the data generation process of the easy TSPTW dataset cannot theoretically guarantee the existence of a feasible solution for each test instance. However, LKH3 find (at least) one feasible solution for each test instance. These results showcase that feasible solutions consistently exist for the synthetic test instances (or datasets) used in this paper, with any reported infeasibility arising solely from the method itself. As mentioned in Section 4, exiting masking mechanism fails to handle problems with complex interdependent constraints (e.g., TSPTW and TSPDL), since it may cause irreversible solution infeasibility during solution construction. Here, we provide a detailed explanation of how our preventative infeasibility masking helps twist the irreversible solution infeasibility, as illustrated in Figure 6 . Given a 5-node (v 0 , v 1 , v 2 , v 3 , v 4 ) TSPTW instance with the time window {[0, 7], [1, 4] , [5, 7] , [2, 5] , [4, 7] }, and the current partial solution v 0 → v 1 , we derive the preventative infeasibility mask by assuming that one of the candidates (v 2 , v 3 , v 4 ) is visited and check the accessibility of assumed arrival time for remaining unvisited nodes. As shown in the left panel, we assume that v 2 is the next visited node and see what would happen if we travel from v 0 → v 1 → v 2 to remaining unvisited nodes v 3 and v 4 . We notice a constraint violation on the assumed tour v 0 → v 1 → v 2 → v 3 since the assumed current time is already 5 at node v 2 ; thus, the earliest arrival time to v 3 is 7, which falls outside the time window of v 3 . Note that another assumed tour v 0 → v 1 → v 2 → v 4 will also violate the time window of v 3 in the end, since it further pushes the assumed arrival time at v 3 . Therefore, v 2In contrast to AM and POMO, the encoder of GFACS outputs edge embeddings, and the decoder consists of a 3-layer multilayer perception (MLP). Each edge embedding h e is processed through:EQUATIONwhere l ∈ [0, 2] denotes the layer index, W l represents a linear layer, and ζ(•) is an activation function, i.e. SiLU for l = 0, 1 and Sigmoid for l = 2. Given the significant correlation between preventative infeasibility masking and the current partial solution τ t , we convert its NAR decoder to an AR one by adding an extra input layer that integrates the embedding of the last edge h et in τ t and the current solution feature f s into the edge embedding itself, i.e., h 0 e = [h e , h et , f s ]. Note that GFACS only consider n 5 (nearest) edges for each node. But for complex constrained problems (e.g., TSPTW), travelling distance may not be sufficient. Instead of using distance as the sole criterion, we select the n 5 neighbors for each node based on the extent of time window overlap, defined as (min(u i , u j ) -max(l i , l j )) for all customer nodes. For the depot, since time window overlap is less relevant, we select n 5 nodes with the earliest l i .We mark the model with Lagrangian multiplier objective function as * (e.g., AM*), and the model further with the preventative infeasibility masking as PIP (e.g., AM*+PIP). We follow their original setups (e.g., model architectures and hyper-parameters) in AM [4] , POMO [5] and GFACS [10] .Pertaining to the PIP-D model (e.g., AM*+PIP-D), we employ an auxiliary decoder (i.e., PIP decoder) that is trained with the ground-truth PIP labels in a supervised manner. To balance the trade-off between training efficiency and empirical performance, we update it periodically. Specifically, within E total training epochs, the PIP decoder is first trained with E init epochs, and then periodically updated E u epochs per E p epochs. To boost the performance, we switch E u to E l for final E l epochs. The detailed settings are presented in Table 9 . For the training epochs that utilize the outputs of the offline (i.e., freezed) PIP decoder as the preventative infeasibility masks, we employ the best-so-far PIP decoder that achieves the highest accuracy on feasible samples, since inaccurate predictions can lead to the exclusion of some feasible candidates. We set the Lagrangian multiplier λ to 1 in the main experiments, with further analyses presented in Appendix D.1. Baselines. We compare our proposed PIP framework against the following baselines, with implementation details provided below:• LKH3 [71], a strong solver designed for a wide range of VRP variants, which we use to generate the (near-)optimal solutions for the test instances with 10,000 trails and 1 run.• OR-Tools [72] , a more flexible solver allowing different combinations of diverse constraints, which we employ the local cheapest insertion as the first solution strategy and the guided local search as the local search strategy with time limit µ for each instance (i.e., 20s for n = 50 and 40s for n = 100 following [27] ). As for TSPDL, we have tried all first solution strategies outlined in the official documentation, yet we still fail to find any feasible solution.• Greedy Heuristics, a classical hand-crafted method considering the local optimal candidates at each step. The Greedy-L heuristic selects the candidate with the shortest distance, and the Greedy-C heuristic selects a node based on the satisfaction of constraints, which is the candidate with the soonest time window end w.r.t current time in TSPTW and the candidate with the minimal draft limit in TSPDL.• AM [4] , a milestone neural AR constructive solver leveraging the Transformer architecture. We implement TSPTW and TSPDL following its default settings.• POMO [5] , an enhanced constructive solver upon AM by considering the symmetry property of VRP solutions. POMO shares a similar architecture with AM. We implement TSPTW and TSPDL following its default settings except removing its stipulated starting node due to the unsuitability to the complex constrained problems, which is also noted in [22] .• GFACS [10] , a neural NAR constructive solver, which introduces the generative flow networks (GFlowNets) to improve the canonical ant colony optimization (ACO) algorithm. It incorporates a heuristic matrix and a pheromone matrix, where the former is parameterized with a neural network, and the latter is updated based on the exploration of multiple ants.We implement TSPTW on it following its default settings.• JAMPR [73] , a state-of-the-art model for a similar problem variant VRPTW, which was further adapted by [9] to solve TSPTW. In this paper, we directly report its result listed in [9] due to the unavailability of source code.• MUSLA [9] , a recent method with designs tailored for TSPTW, incorporating problemspecific features and a large supervised learning dataset, where OSLA is its one-step version and MUSLA adapt adopts an adaptive inference strategy to balance the optimality gap and the solution feasibility. Results are adopted from [9] due to the unavailability of source code, with our 'Easy' setting corresponding to their 'Medium' one in [9] . Our PIP leverages the Lagrangian multiplier method to guide neural policy search. While some existing methods, such as [8] , employ bi-level optimization techniques to jointly update the Lagrangian multiplier and the primal objective function, we find these approaches to be inefficient. Therefore, we opt to fix the value of the Lagrangian multiplier λ and focus solely on optimizing the primal variables in this paper. Here, we investigate the effect of different Lagrangian multipliers on the solution quality and feasibility. The results are shown in the left panel of Figure 7 , where we observe that a larger λ results in a lower infeasible rate but with an increased optimality gap, whereas a smaller λ reduces the optimality gap but at the expense of a much higher infeasible rate. Consequently, we set λ = 1 for a balance, prioritizing the optimization of the primal objective function.We report the average objective values and optimality gaps in Tables 1 and 2 . It should be noted, however, that these metrics are calculated across different sets of feasible instances. Considering that a complete overlap of feasible instances across all baselines is impractical (due to the 100% infeasibile rates of some baselines), we present additional results on a set of overlapped feasible instances across POMO variants to facilitate a more conprehensive comparison. As shown in Table 10 , the models with our PIP consistently outperforms POMO in terms of solution quality. On the easy dataset, PIP-D exceeds PIP due to the incorporation of data augmentation during inference. In summary, despite the different instance sets, the optimality gaps displayed in Table 1 (w/o overlap) exhibit numerical patterns and conclusions that align with those empirically observed on the overlapped sets. As illustrated in Figure 1 , the feasibility masking (i.e., n-step PIP) in complex constraints is NP-Hard.While iterating over all future possibilities would make PI masking complete, it is computationally inefficient. Therefore, we approximate it with one-step PI masking, whose efficiency is validated in Table 4 and Table 5 . To enhance training efficiency, we use an auxiliary decoder to further approximate one-step PI masking, avoiding the need to acquire it continuously during training.Besides, we further explore some strategies for accelerating the training and inference of PIP.• Apply sparse strategies to refine PIP calculations. Due to the O(n 2 ) complexity of PIP, applying it to all the unvisited candidate nodes will be computationally expensive. For large-scale problems, we only consider top K neighbours, which is implemented on GFACS.Results in Table 3 show that GFACS*+PIP-D maintains similar training and inference times as GFACS* on TSPTW-500 (i.e., 28.3h and 6.5m vs. 28.1h and 6.4m).• Couple with the state-of-the-art solvers (e.g. LKH3). Our PIP is empirically verified to be efficient due to its capability to obtain good and feasible solutions within a very short time (LKH3: 1.4d vs POMO*+PIP-D: 48s), while LKH3 can get near-optimal solutions with prolonged time. To leverage the strengths of both approaches, we use our PIP-D to provide better initial solutions for LKH3. As shown in Table 11 , this combination reduces the infeasibility rate from 53.11% to 0.21% and improves the objective from 51.65 to 51.25 within only a few seconds per instance. Notably, initializing LKH3 with POMO*+PIP-D outperforms the default LKH3 setup (10,000 trials), achieving slightly better solution quality while using only 27% of the inference time (9 hours vs. 1.4 days). We also show the progress of objective value and instance-level infeasibility rate over inference times in Figures 8, 9 for clearer comparison.• Fine-tune Lagrangian method (*) with PI masking. On top of the basic Lagrangian method (e.g., POMO*), PIP further employs the preventative infeasibility (PI) masking throughout the training process. We would like to note that there is another way to exploit the preventative infeasibility information and reduce the computational complexity, i.e., by leveraging the PI mask to fine-tune the pretrained Lagrangian method. The comprehensive empirical results on Medium TSPTW-50 are shown in Table 12 , where the first two methods denote the PIP and Lagrangian methods applied to POMO, respectively. The results indicate that a few steps of fine-tuning the Lagrangian method using PIP masks can yield favorable improvements in the feasible rate and the optimality gap, and significantly reduce the training complexity as well.• Early stop of PI masking. At each solution construction (i.e., decoding) step, our PIP leverages preventative infeasibility masking to proactively steer the policy search to (near-)feasible regions, leading to increased computational overheads as the problem sizes scale up. Here, we explore the potential of early stopping PIP, where PIP is only employed during the initial steps of solution construction. Based on our empirical observation depicted in the left panel of Figure 10 , infeasibility predominately occurs in the first few steps of the entire process. This observation reveals the possibility of merely acquiring PI masks for the initial few steps, which could improve the training efficiency. We leave it to future work. We use an auxiliary decoder (i.e., PIP decoder) in the proposed PIP-D framework, whose goal is to learn and predict PI masks by identifying infeasible candidates based on the current partial solution.We formulate it as a binary classification task. Here, we demonstrate the efficacy of the learned PIP decoder through various evaluation metrics, including accuracy, recall (of infeasible samples), and specificity (recall of feasible samples). As shown in the last two panels of Figure 10 , our PIP decoder can accurately predict the PI masks across all hardness levels in TSPTW and TSPDL, especially for the more complex constrained one. This indicates that our PIP decoder aligns well with the goal of ensuring high accuracy in feasibility predictions.D.5 Performance under different inference time budget.For a fair comparison, we also extend the inference time (by sampling more solutions and data augmentation) of the baselines to a similar one as POMO + PIP and POMO + PIP-D. Results show that incorporating the Lagrangian multiplier (POMO* and AM*) may lead to some improvement (in Table 1 and 2 ), but not for the cases in Table 13 under complex constraints and larger scales. Even with prolonged inference time, existing methods do not deliver any feasible solutions for the studied complex constrained VRP. In contrast, our PIP-D with 48s time significantly reduces infeasibility from 100% to 6.48% compared to baselines running for 2.5m, and exhibits an optimality gap of around only 0.3%. Furthermore, PIP can perform even better with more inference time. 11 , the performance disparity among different methods is significant across all hardness levels, especially in more complex constrained problems, underscoring the effectiveness of our method. Note that this result further supports the improvement of our PIP framework reported in Table 1 . Performance under different hyper-parameters. We explore the variance in model performance under different hyper-parameters. We fix the total update epoch of the PIP decoder and adjust its interval. Results in This paper focuses on real-world scenarios and proposes a novel Proactive Infeasibility Prevention (PIP) framework to enhance the capabilities of neural methods towards solving more complex VRPs. Potential positive societal impacts include: 1) enhancing industrial efficiency, e.g., in logistics and transportation. By preemptively identifying infeasible solutions, it can reduce computational overheads and improve the efficiency of decision-making process; 2) advancing the AI and operation research (OR) communities. Our PIP framework aims to alleviate the existing challenges in the neural VRP solvers, thereby promoting the advancement of AI as well as OR. On the other hand, negative societal impacts may include environmental unfriendliness due to computational resource usage.The used assets in this work are listed in Table 16 , which are all open-source for academic research. We will release our source code with the MIT License.Table 16 : Used assets, licenses, and their usage.Code LKH3 [71] Available for academic use Evaluation OR-Tools [72] Apache-2.0 license Evaluation AM [4] MIT License Revision POMO [5] MIT License Revision GFACS [10] MIT Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 


Decision: This paper introduces and analyzes a  Proactive Infeasibility Prevention (PIP) framework increasing the number of feasible solutions generated when solving complex VRPs, ultimately leading to better quality final solutions found. The paper has several strengths - the problem addressed is an important one, the PIP approach is novel, and a comprehensive experimental analysis (including the generation of several sets of new results during the rebuttal discussion period) demonstrate the advantage of the approach in dealing with VRP/TSP problems with complex constraints. Finally the paper is well-written and understandable.

There were two main perceived weaknesses of the paper. First, was the fact that the proposed method incurs significantly more computational expense  than previous approaches to addressing this problem. Given this fact, there were concerns that some of the experimental analyses that were performed did not adequately make a fair comparison against prior existing approaches.   These concerns led to the addition of several new sets of experiments during the author rebuttal period, which clarified and further highlighted PIP's ability to outperform existing approaches in solution quality, and also showed the ability to further boost solution quality when used in combination. A second observed weakness was the lack of any theoretical model justifying the PIP approach. Subsequent discussion here pointed out that component elements of the approach (e.g., use of Lagrangian relaxation) did have theoretical justification, and consensus was reached that development of a theoretical basis was an important area for future research. Please address all original reviewer concerns and  incorporate all of these new results in the final version of the paper.

Reviews: review: paper_summary: paper_summary: This paper proposes a Proactive Infeasibility Prevention (PIP) framework to enhance the ability of neural methods to handle complex constraints in Vehicle Routing Problem (VRP).
The PIP framework integrates the Lagrangian multiplier to enhance constraint awareness and introduces preventative infeasibility masking to proactively guide the solution construction process. Additionally, an extended version called PIP-D employs an auxiliary decoder and two adaptive strategies to learn and predict masking information, reducing computational costs during training.
These methods were extensively tested on different levels of constraint hardness in the Traveling Salesman Problem with Time Window (TSPTW) and Traveling Salesman Problem with Draft Limit (TSPDL) variants. The results demonstrate that the proposed methods enhance the capabilities of neural methods, significantly reducing infeasibility rates and improving solution quality.
strength_weakness: strength_weakness: 1. Generality of the Trained Model: It is unclear whether the trained model generalizes well to unseen instances, raising concerns about its robustness and applicability to different problem settings.
2. Lack of Theoretical Justification for PIP: The paper does not provide a theoretical justification for the superiority of using the Proactive Infeasibility Prevention (PIP) method, leaving its theoretical advantages unproven.
3. Lack of Classical methods with time limit: There are no comparison with the classical methods (e.g., LKH3) with a time limit. I am concerned that LKH3 could find "good" feasible solutions within a short time (e.g., within 5 minutes).
questions: questions: - Can the proposed idea be applied to methods other than constructive methods? If it can be used to guide the search to avoid infeasibility, it would be a more practical idea.
- Alternatively, after the (possibly infeasible) solution constructed by a constructive method, is it easy to recover a feasibility through neighborhood search? If so, can you compare this with the avoidance of infeasibility using PIP?
- Is the model constructed with PIP tailored to the training dataset? In other words, when using the model constructed with PIP for instances in different area, does the number of infeasible instances decrease compared to the model without PIP?
- Is it possible to use the PIP framework for more general VRP (e.g., with multiple vehicles or capacity constraints)?
- Are there any experimental comparisons between the proposed model and LKH3 with time limit?
limitations: limitations: The limitations of this work are discussed in the last section.
paper_summary: paper_summary: The paper  propose a novel Proactive Infeasibility Prevention (PIP) framework to advance the capabilities of neural methods towards more complex VRPs, and further investigates the Lagrange multiplier method for soft objective in VRPs, presenting the PIP (& PIP-D) for the hard case where Lagrange multiplier method difficult to find feasible solution. And the experiments show that the PIP is prior to Lagrange multiplier method and has generality to TSPTW and TSPDL.
strength_weakness: strength_weakness: The author's method actually addresses the shortcomings of the Lagrange multiplier method, but it is not mentioned in the Introduction, and the problems raised in the Introduction can actually be solved only using the Lagrange multiplier method.   
The method involves complex calculations, particularly when integrating the Lagrangian multiplier and PI masking, which can be computationally intensive​. 
The paper does not thoroughly address the scalability of the proposed method when applied to extremely large datasets, which may limit its practical application in some real-world scenarios​​.
questions: questions: [1]What are the potential challenges and limitations in scaling the PIP method to larger and more complex combinatorial optimization problems?

[2]How can we further optimize or modify the proposed method to reduce computational complexity without sacrificing solution quality or feasibility?
limitations: limitations: yes
paper_summary: paper_summary: This paper addresses the challenge of predicting feasible solutions for VRPs with complex constraints. It introduces two novel methods to enhance existing algorithms: i) integrating constraints directly into the optimization objective using the Lagrangian Multiplier Approach; ii) PIP framework, which proactively excludes potential nodes during solution prediction to prevent future infeasibility. Besides, a neural decoder PIP-D is proposed to reduce computational complexities.
Experiments on two VRP variants demonstrates significant advancements in feasible prediction compared to the baselines.
strength_weakness: strength_weakness: Major concerns:

- **Computational cost.** From the numerical results, PIP incurs approximately twice the computational time compared to baseline neural solvers (e.g., AM and POMO). Despite attempts to mitigate this by introducing neural decoder PIP-D, the reported results don't exhibit a corresponding reduction in time. The authors may clarify this more.

- **Lack of comparisons** between different step numbers in PIP. The authors only conducted experiments of one-step PIP. It is crucial to include comparisons with different step numbers to assess the impact of the "looking ahead" masking strategy. For instance, additional results from zero-step PIP (directly masking nodes violating constraints) and two-step PIP. Here, I emphasize zero-step PIP since its masking strategy differs from the ones of the baseline neural solvers, as I understand it. 

- Only two TSP datasets are considered. The paper is titled "Learning to Handle Complex Constraints for **Vehicle Routing Problems (VRPs)**", however both datasets are **variants of TSPs**. The key difference between VRPs and TSPs is **capacity constraints**. I would like to see VRPs with various complex constraints (such as time windows, pick-up and delivery, split delivery etc.), otherwise, the authors should change the title and consider more TSP variants. 

- Heuristic algorithms for TSPs and VRPs with complex constraints have been well developed (such as LKH3). I do not see that the proposed algorithm outperforms the state-of-the-art heuristic algorithms (not neural algorithms).
questions: questions: - What's the time cost of PIP? Is it $time(POMO^*+PIP)$ - $time(POMO^*)$?

- What are the masking strategies of baselines AM and POMO? Do they only mask visited nodes, ignoring other constraints?

- What is the stopping criteria for each algorithm? I know that running LKH3 for 1 second would often generate almost optimal solutions to TSPs with 100 customers, hence there is no need to run it for 14 hours, as shown in Table 2. I recommend the authors use the same time limit (and usually a few seconds for small-sized TSPs considered in this paper) for each baseline algorithm and even present the progress of the objective value for the best solutions found by each algorithm. This could ensure a fair comparison among them.
limitations: limitations: The author has discussed the limitations in section 6: one potential limitation is that it may not improve performance on all backbone solvers and all VRP variants.
paper_summary: paper_summary: This paper introduces a novel framework called Proactive Infeasibility Prevention (PIP) to enhance the capability of neural methods in addressing complex Vehicle Routing Problems (VRPs) with interdependent constraints, such as the Traveling Salesman Problem with Time Window (TSPTW) and TSP with Draft Limit (TSPDL). The authors propose integrating the Lagrangian multiplier method and preventative infeasibility masking into the solution construction process to improve both the feasibility and optimality of solutions. Additionally, an enhanced version, PIP-D, employs an auxiliary decoder to predict infeasibility masks, potentially reducing computational costs during training. The paper presents extensive experiments demonstrating the effectiveness of PIP in reducing infeasible rates and improving solution quality.
strength_weakness: strength_weakness: The ideas of using the Lagrangian multiplier [1] and preventive mask functions[2] are not new.

The code is not provided, which limits the reproducibility of the work.

It's not clear how to obtain the training label of infeasibility mask.  The computation could be expensive to get a global and accurate infeasibility mask.

Some related work from recent AI and OR fields are missing.

The paper could benefit from additional experiments on a wider range of VRP variants and real-world datasets to further validate the generalizability of the proposed methods.

The presentation could be improved by providing more context on how this work relates to and advances the current state-of-the-art in neural VRP solvers.

The computational efficiency of the PIP-D framework should be more thoroughly analyzed, especially in terms of how it scales with problem size.

The paper may lack a deeper discussion on the theoretical underpinnings of the proposed methods and their potential implications for the broader field of combinatorial optimization.

[1] Qiaoyue Tang, Yangzhe Kong, Lemeng Pan, and Choonmeng Lee. Learning to solve soft-constrained vehicle routing problems with lagrangian relaxation. arXiv preprint arXiv:2207.09860, 2022.

[2] Hou, Qingchun, et al. Generalize learned heuristics to solve large-scale vehicle routing problems in real-time. The Eleventh International Conference on Learning Representations. 2023.
questions: questions: How does the PIP framework generalize to other VRP variants beyond TSPTW and TSPDL?

Can the authors provide more insights into the computational efficiency of PIP-D, especially regarding its scalability?
limitations: limitations: The authors have acknowledged some limitations, such as the potential inability of PIP to improve performance on all backbone solvers and VRP variants. However, the paper could benefit from a more detailed discussion on the robustness of the PIP framework when faced with different levels of constraint hardness in real-world scenarios.

Comments: title: Thank you for your support!
Dear Reviewers DEtT and Seke:
Thank you for your support!
Official Comment by Reviewer azSu
Thank you for your support and we are happy to have addressed your concerns!
Further results for comprehensive comparison (4/4)
Further results for comprehensive comparison (3/4)
Further results for comprehensive comparison (2/4)
Further results for comprehensive comparison (1/4)
Thanks for the comments and please kindly check our further clarifications (continued)
Thanks for the comments and please kindly check our further clarifications.
comment: Dear Reviewer Seke,
####
The discussion period is approaching its end (in less than 12 hours). Please kindly let us know if our response resolves your concerns. Regarding your concern of our code release, we have also forwarded our code to the Area Chair. We would greatly appreciate it if you could give us any feedback. 
####
Thank you again for your valuable comments and suggestions.
####
Best Regards,

Authors
We greatly appreciate your time in reviewing our paper and reading the follow-up rebuttal! We're thrilled for your recognition of our work, and will include the additional analysis in the revised paper. Thanks a lot!
Thank you very much for your thoughtful response and the additional analysis. Most of my concerns have been resolved, and my score has been raised.
The authors have provided extensive comments and new results in response to the criticisms raised in your reviews. Has this response addressed your main concerns? If not, are there additional questions you would like to pose? The author discussion period is scheduled to end tomorrow (Aug 15). Please respond.
We really appreciate your constructive comments and strong support for our paper. Thank you!
Thank you very much for your thoughtful response and the additional analysis. Most of my concerns have been resolved.

I understand that the theoretical analysis of PIP and its application to iterative methods are considered future work.
I agree that the proposed methods and its experiments presented in this paper are valuable, and I raised my score.
Dear Reviewer Seke,

Thank you once again for your insightful comments and helpful suggestions. As the deadline for author-reviewer discussions is approaching, we would greatly appreciate it if you could take a moment to review our rebuttal. Please let us know if you have any further questions or concerns. Thank you very much for your time.

Best,
Authors
Dear Reviewer DEtT,

Thank you once again for your insightful comments and helpful suggestions. As the deadline for author-reviewer discussions is approaching, we would greatly appreciate it if you could take a moment to review our rebuttal. Please let us know if you have any further questions or concerns. Thank you very much for your time.

Best,  
Authors
Dear Reviewer azSu,

Thank you once again for your insightful comments and helpful suggestions. As the deadline for author-reviewer discussions is approaching, we would greatly appreciate it if you could take a moment to review our rebuttal. Please let us know if you have any further questions or concerns. Thank you very much for your time.

Best,  
Authors
We really appreciate your insightful suggestions and we will include the additional results in the revised paper. Thank you!
Thanks for addressing my comments. Now I'm satisfied with the comparison and raise my score. I expect the authors to include those results in the updated manuscript and edit part of their conclusions (the proposed approach vs LKH).
**Table R8**: Results on POMO*+PIP(-D)+LKH3 with the similar $\underline{\text {total inference time}}$ limit as doubled time used in POMO*+PIP(-D).

|              |                  |           |                                 |    $n=50$     |                  |              |                                 | $n=100$   |                  |   
|:------------:|:----------------:|:---------:|:-------------------------------:|:-------------:|:----------------:|:------------:|:-------------------------------:|:---------:|:----------------:|   
| **Hardness** |    **Method**    |           | $\underline{\text{Total Time}}$ |   **Obj.**    | **Inst. Infsb%** |              | $\underline{\text{Total Time}}$ | **Obj.**  | **Inst. Infsb%** |   
|     Easy     |       LKH3       |           |               46s               |     7.57      |      0.17%       |              |              1.8m               |   11.17   |      2.54%       |   
|     Easy     |  POMO*+PIP+LKH3  |           |               48s               |     7.48      |      0.20%       |              |              1.8m               |   10.46   |      0.01%       |   
|     Easy     | POMO*+PIP-D+LKH3 |           |               48s               |     7.47      |      0.32%       |              |              1.8m               |   10.52   |      0.00%       |   
|              |                  |           |                                 |               |                  |              |                                 |           |                  |   
|    Medium    |       LKH3       |           |               46s               |     13.36     |      0.67%       |              |              1.8m               |   20.42   |      26.66%      |   
|    Medium    |  POMO*+PIP+LKH3  |           |               48s               |     13.37     |      1.12%       |              |              1.9m               |   19.38   |      0.20%       |   
|    Medium    | POMO*+PIP-D+LKH3 |           |               48s               |     13.41     |      0.77%       |              |              1.9m               |   19.50   |      0.05%       |   
|              |                  |           |                                 |               |                  |              |                                 |           |                  |   
|     Hard     |       LKH3       |           |               44s               |     23.87     |      99.09%      |              |              1.9m               |     /     |     100.00%      |   
|     Hard     |  POMO*+PIP+LKH3  |           |               44s               |     25.65     |      2.99%       |              |              1.8m               |   51.39   |      15.80%      |   
|     Hard     | POMO*+PIP-D+LKH3 |           |               44s               |     25.68     |      3.23%       |              |              1.8m               |   51.37   |      6.24%       |   

####
####

**Table R9**: Results on LKH3 with the similar $\underline{\text {instance inference time}}$ limit as GFACS*+PIP(-D) (+LKH3) on $n=500$ (first 10 instances).


|                   |                |     $n=500$      |              |
| :---------------: | :------------: | :--------------: | :----------: |
|    **Method**     | **Inst. Time** | **Inst. Infsb%** | **Gap to #** |
|  LKH3 (Default)   |      26m       |      0.00%       |      #       |
|       LKH3        |      6.5m      |      0.00%       |    0.59%     |
|      GFACS*       |      6.4m      |      57.81%      |    21.32%    |
|   GFACS*+PIP-D    |      6.5m      |      0.00%       |    11.52%    |
|       LKH3        |      13m       |      0.00%       |    0.25%     |
| GFACS*+PIP-D+LKH3 |      13m       |      0.00%       |    1.10%     |
**Table R7**: Results on POMO*+PIP(-D)+LKH3 with the similar $\underline{\text {instance inference time}}$ limit as doubled time used in POMO*+PIP(-D).
                                                                                                                                                                    
|              |                  |                                 |  $n=50$  |                  |                                 | $n=100$  |                  | 
|:------------:|:----------------:|:-------------------------------:|:--------:|:----------------:|:-------------------------------:|:--------:|:----------------:| 
| **Hardness** |    **Method**    | $\underline{\text{Inst. Time}}$ | **Obj.** | **Inst. Infsb%** | $\underline{\text{Inst. Time}}$ | **Obj.** | **Inst. Infsb%** | 
|     Easy     |       LKH3       |              0.77s              |   7.33   |      0.00%       |              1.8s               |  10.31   |      0.00%       | 
|     Easy     |  POMO*+PIP+LKH3  |              0.77s              |   7.33   |      0.00%       |              1.8s               |  10.29   |      0.00%       | 
|     Easy     | POMO*+PIP-D+LKH3 |              0.77s              |   7.33   |      0.00%       |              1.8s               |  10.31   |      0.00%       | 
|              |                  |                                 |          |                  |                                 |          |                  | 
|    Medium    |       LKH3       |              0.77s              |  13.04   |      0.00%       |              1.8s               |  18.89   |      0.00%       | 
|    Medium    |  POMO*+PIP+LKH3  |              0.77s              |  13.05   |      0.00%       |              1.8s               |  18.91   |      0.00%       | 
|    Medium    | POMO*+PIP-D+LKH3 |              0.77s              |  13.05   |      0.00%       |              1.8s               |  18.92   |      0.00%       | 
|              |                  |                                 |          |                  |                                 |          |                  | 
|     Hard     |       LKH3       |              0.76s              |  25.56  |      8.81%       |              1.8s               |  50.66   |     75.54%       | 
|     Hard     |  POMO*+PIP+LKH3  |              0.76s              |  25.61   |      0.12%       |              1.8s               |  51.24   |      1.63%       | 
|     Hard     | POMO*+PIP-D+LKH3 |              0.76s              |  25.61   |      0.05%       |              1.8s               |  51.25   |      0.36%       |
**Table R6**: Results on LKH3 with the similar $\underline{\text {total inference time}}$ limit as POMO*+PIP(-D).

                                                                                                                                                                 
|              |                |                                 |  $n=50$  |                  |                                 | $n=100$  |                  |
|:------------:|:--------------:|:-------------------------------:|:--------:|:----------------:|:-------------------------------:|:--------:|:----------------:|
| **Hardness** |   **Method**   | $\underline{\text{Total Time}}$ | **Obj.** | **Inst. Infsb%** | $\underline{\text{Total Time}}$ | **Obj.** | **Inst. Infsb%** |
|     Easy     | LKH3 (Default) |              4.6h               |   7.31   |      0.00%       |              8.5h               |  10.21   |      0.00%       |
|     Easy     |      LKH3      |               26s               |   8.81   |      99.29%      |               58s               |    /     |     100.00%      |
|     Easy     |   POMO*+PIP    |               21s               |   7.50   |      0.00%       |               48s               |  10.57   |      0.00%       |
|     Easy     |  POMO*+PIP-D   |               21s               |   7.49   |      0.00%       |               48s               |  10.66   |      0.00%       |
|              |                |                                 |          |                  |                                 |          |                  |
|    Medium    | LKH3 (Default) |               7h                |  13.02   |      0.00%       |              10.8h              |  18.74   |      0.00%       |
|    Medium    |      LKH3      |               25s               |  13.05   |      39.91%      |               63s               |    /     |     100.00%      |
|    Medium    |   POMO*+PIP    |               21s               |  13.40   |      0.90%       |               48s               |  19.61   |      0.19%       |
|    Medium    |  POMO*+PIP-D   |               21s               |  13.45   |      0.65%       |               48s               |  19.79   |      0.03%       |
|              |                |                                 |          |                  |                                 |          |                  |
|     Hard     | LKH3 (Default) |               7h                |  25.61   |      0.12%       |              1.4d               |  51.24   |      0.07%       |
|     Hard     |      LKH3      |               22s               |    /     |     100.00%      |               54s               |    /     |     100.00%      |
|     Hard     |   POMO*+PIP    |               21s               |  25.66   |      2.67%       |               48s               |  51.42   |      16.27%      |
|     Hard     |  POMO*+PIP-D   |               21s               |  25.69   |      3.07%       |               48s               |  51.39   |      6.48%       |
**Table R5**: Results on LKH3 with the similar $\underline{\text {instance inference time}}$ limit as POMO*+PIP(-D).

  
|              |                |                                 |  $n=50$  |                  |                                 | $n=100$  |                  | 
|:------------:|:--------------:|:-------------------------------:|:--------:|:----------------:|:-------------------------------:|:--------:|:----------------:| 
| **Hardness** |   **Method**   | $\underline{\text{Inst. Time}}$ | **Obj.** | **Inst. Infsb%** | $\underline{\text{Inst. Time}}$ | **Obj.** | **Inst. Infsb%** |  
|     Easy     | LKH3 (Default) |               27s               |   7.31   |      0.00%       |               49s               |  10.21   |      0.00%       | 
|     Easy     |     LKH3       |              0.37s              |   7.35   |      0.00%       |              0.9s               |  10.37   |      0.00%       | 
|     Easy     |   POMO*+PIP    |              0.38s              |   7.50   |      0.00%       |              0.9s               |  10.57   |      0.00%       | 
|     Easy     |  POMO*+PIP-D   |              0.38s              |   7.49   |      0.00%       |              0.9s               |  10.66   |      0.00%       | 
|              |                |                                 |          |                  |                                 |          |                  | 
|    Medium    | LKH3 (Default) |               40s               |  13.02   |      0.00%       |              1.0m               |  18.74   |      0.00%       | 
|    Medium    |      LKH3      |              0.37s              |  13.06   |      0.00%       |              0.9s               |  19.00   |      0.00%       | 
|    Medium    |   POMO*+PIP    |              0.38s              |  13.40   |      0.90%       |              0.9s               |  19.61   |      0.19%       | 
|    Medium    |  POMO*+PIP-D   |              0.38s              |  13.45   |      0.65%       |              0.9s               |  19.79   |      0.03%       | 
|              |                |                                 |          |                  |                                 |          |                  | 
|     Hard     | LKH3 (Default) |               40s               |  25.61   |      0.12%       |              3.2m               |  51.24   |      0.07%       | 
|     Hard     |      LKH3      |              0.37s              |  25.43   |      30.60%      |             0.9s                |  49.94   |       97.28%     | 
|     Hard     |   POMO*+PIP    |              0.38s              |  25.66   |      2.67%       |              0.9s               |   51.42  |      16.27%      | 
|     Hard     |  POMO*+PIP-D   |              0.38s              |  25.69   |      3.07%       |              0.9s               |  51.39   |      6.48%       |
Additionally, we would like to clarify that:

1. **The motivation of this work is to highlight the issue of “NP-hard masking” in neural solvers and enhance neural solvers with the ability to handle complex constraints**, as comprehensively demonstrated in our main paper (Tables 1, 2, and 3). Thus, outperforming LKH3 is not our primary objective or key contribution. Instead, our proposed framework has significantly reduced the infeasibility rate and substantially improved solution quality for both autoregressive and non-autoregressive neural solvers. These advancements establish our framework as the new state-of-the-art in the domain of neural solvers.

2. **Most existing NCO solvers [2-8] use LKH3 with the default settings from [1]** (including maximum trials, number of runs, and parallelism settings) to estimate the optimality gap and provide a runtime reference. The primary focus is on understanding the performance gap between neural solvers [10, 11] and traditional solvers, rather than strictly outperforming it.

3. We acknowledge that **neural solvers have yet to fully replace traditional solvers like LKH3, but rather catching up with it step by step**. We believe that we provide insights for complex constraints handling for NCO, which is a step towards practical applications. Also, we would like to emphasize that our PIP(-D) is a **generic framework**, making our research orthogonal to the development of new, powerful neural solvers that have the potential to fully outperform LKH in future works. 

We appreciate the reviewer for helping our work to be more rigorous. The new results will be incorporated into the revised manuscript, along with a deeper discussion of these points. Due to the rebuttal limit, we will plot the progress curves (like Figure R2 and R3 in the attached PDF) based on the results of all the new tables for a more intuitive illustration. Also, **our code has already been forwarded to the AC**. We hope our discussions are now comprehensive. Please let us know if you have any remaining concerns.

####
####

```
References:
[1] Attention, learn to solve routing problems!. ICLR, 2018.
[2] Reinforcement learning for solving the vehicle routing problem. NeurIPS, 2018.
[3] POMO: Policy optimization with multiple optima for reinforcement learning. NeurIPS, 2022.
[4] Simulation-guided beam search for neural combinatorial optimization. NeurIPS, 2022.
[5] Learning to search feasible and infeasible regions of routing problems with flexible neural k-opt. NeurIPS, 2023.
[6] BQ-NCO: Bisimulation quotienting for generalizable neural combinatorial optimization. NeurIPS, 2023.
[7] Neural combinatorial optimization with heavy decoder: Toward large scale generalization. NeurIPS, 2023.
[8] MVMoE: Multi-Task Vehicle Routing Solver with Mixture-of-Experts. ICLM, 2024.
[9] Machine learning for combinatorial optimization: a methodological tour d’horizon. European Journal of Operational Research, 2021.
[10] Graph-based diffusion solvers for combinatorial optimization. NeurIPS, 2023.
[11] MOCO: A Learnable Meta Optimizer for Combinatorial Optimization. arXiv preprint, 2024.
```

#### 

Note for Tables R7 & R8: For a fair comparison, we kept the POMO*+PIP(-D) models (with inference time denoted as $T_{PIP}$) unchanged and allowed LKH3 to search for a similar duration as $T_{PIP}$ to showcase the performance of POMO*+PIP(-D)+LKH3 (i.e., the total inference time is $2 \times T_{PIP}$). LKH3 was also given the same time limit for consistency.
Thank you for the discussion and insightful comments!

We are pleased to have addressed most of your concerns. Regarding the remaining concern about a comprehensive comparison between our PIP and LKH3, we have provided results on Hard TSPTW-100 in the previous rebuttal given the constraints of the rebuttal period and PDF page limit. Nevertheless, we understand that this may not be comprehensive enough. Following your suggestion, we now conducted additional experiments using LKH3 with the same time limits as our proposed approach **across various instance hardness levels (Easy/Medium/Hard) and scales (50/100/500 nodes)**. The time limits were set in two ways for comprehensive evaluation: similar $\underline{\text {instance inference time}}$ without parallelization (new Tables R5, R7, and R9) and similar $\underline{\text {total inference time}}$ with parallelization (new Tables R6 and R8). We hope these new empirical results offer a more comprehensive picture.

**[A. For $n=50$ and $100$]**  Firstly, we follow your suggestion and only allow LKH3 to run a similar inference time as our approach $\underline{\text {per instance}}$. As shown in **Table R5**, POMO*+PIP(-D) outperforms LKH3 on Hard datasets, while maintaining competitive results on Easy and Medium datasets. While comparing per-instance time might seem fair for CPU-based LKH3, ignoring parallelization could disadvantage GPU-based solvers (thus not fair for GPU-based solvers). To further explore this, we conducted another experiment but with the similar $\underline{\text {total inference time}}$ limit across both methods. Results, in **Table R6**, show that our POMO*+PIP(-D) performs consistently better than LKH3 across most of the hardness. We recognize the challenge of achieving absolute fairness when comparing CPU-based and GPU-based solvers, and we note that total time comparison is a common practice in most existing NCO papers (e.g. those in [1-8]). We appreciate the reviewer for providing such an insightful comment, which highlights the need for developing new metrics that facilitate fairer comparisons between neural and traditional solvers—a promising direction for future work.

Moreover, we agree with the reviewer that LKH3 is indeed a robust and well-developed solver. However, neural solvers offer unique advantages, such as highly efficient parallelization and reduced reliance on hand-crafted heuristic rules [9]. To leverage the strengths of both approaches and further reveal the practical usage of our method, we explored a hybrid method that combines our PIP(-D) framework with LKH3. Similarly, we provide the comprehensive experimental results conducted across various instance hardness levels (Easy/Medium/Hard) and scales (50/100/500 nodes). The results, as shown in **Table R7 and R8**, reveal that **LKH3’s search efficiency can be significantly enhanced when initialized with solutions from our PIP(-D) framework**.

**[B. For $n=500$]** As shown in **Table R9**, we acknowledge that our PIP(-D) implemented on GFACS* currently underperforms LKH3. However, the integration of GFACS*+PIP(-D) with LKH3 has already made significant strides in closing this performance gap, particularly in enhancing feasibility handling compared to using the GFACS backbone alone. We would like to clarify that the primary motivation for these experiments here is to demonstrate the scalability and broad applicability of our PIP framework in enhancing various backbone neural models. While we acknowledge that PIP may not elevate every model to state-of-the-art performance, which is both reasonable due to _“No Free Lunch”_, our goal is to show its potential to significantly boost performance across a wide range of neural models.
Thanks for addressing my comments. I still think the comparison between the proposed approach and LKH3 is not comprehensive enough (e.g., easy/medium/hard instances with more than 100 nodes, and the 500-customer instances in Section 5.2; **with the same time limit as the proposed approach**) since **LKH3 is the most important baseline algorithm** in my opinion. For example, from Figure R3, we may conclude that, for those hard instances, the proposed approach could obtain high-quality solutions quite fast but LKH3 could obtain solutions of better quality with a relatively longer time. But from Section 5.2, if we run LKH3 for a few minutes (say 6.5 min), we may already achieve a very small gap (e.g., smaller than 11.95%). That is why I said the comparison in the previous manuscript is not fair.