Efficient Neural Neighborhood Search for Pickup and Delivery Problems

1 Introduction

Efficient neighborhood search functions as the key component of powerful heuristics for pickup and delivery problems (PDPs) [Parragh et al., 2008]. It involves an iterative search process which transforms a solution into another candidate in its current neighborhood, hopefully in an efficient way. Designing the neighborhood and search rules usually determines the efficiency and even the success of a solver. However, they are often problem-specific and manually engineered with a lot of trial and error, which needs to redesign when changes occur in constraints or objectives. These limitations may hinder the applications in the rapidly evolving industries.

On the other hand, recent neural methods for vehicle routing problems (VRPs) have emerged as promising alternatives to traditional heuristics (e.g., [Ma et al., 2021b]). They are usually faster, and more importantly, could automate the design of heuristics for new variants where no hand-crafted rule is available [Kwon et al., 2020]. However, the prevailing neural methods mainly focus on travelling salesman problem (TSP) or capacitated vehicle routing problem (CVRP), where efficient solvers for PDPs are rarely studied. The PDP is ubiquitous in logistics, robotics, meal-delivery services, etc [Parragh et al., 2008], which optimizes the route for pickupdelivery requests and is characterised by the precedence constraint (pickup before delivery). Despite the first attempt in [Li et al., 2021b], which learns a construction method to build a PDP solution (route) in seconds, it leaves a considerable gap to traditional heuristics in solution quality.

∗Equally contributed.

†Zhiguang Cao is the corresponding author.

1 Code is available at https://github.com/yining043/PDP-N2S. Please refer to https://arxiv.org/abs/2204.11399 for an extended version of this paper with full results of Table 7 and more discussion.

To reduce the gap, we propose an efficient Neural Neighborhood Search (N2S) approach for PDPs, based on a novel Transformer styled policy network with the encoder-decoder structure. To our knowledge, the most similar existing policy network to ours is the DACT in [Ma et al., 2021b]. Also as an improvement method, DACT learns to encode the current solution and transform it into another one using the 2-opt, insert, or swap decoder. Among them, the 2-opt decoder, which considers reversing a segment of the solution, performed the best for TSP and CVRP. However, it is not suitable for PDPs since the precedence constraint can be easily violated by segment inversion. Though the insert decoder performed well on small-scale PDPs, which considers removing and reinserting a single node, its performance significantly drops on largerscale or highly constrained PDPs (see Section 5). Conversely, our N2S tackles the precedence constraint more efficiently by allowing a pair of pickup-delivery nodes to be simultaneously operated in the neighborhood search through two customized removal and reinsertion decoders as shown in Figure 1.

Another challenge lies in the design of the encoders. It was well revealed in [Ma et al., 2021b] that the vanilla Transformer encoder [Vaswani et al., 2017] failed to correctly encode route solutions since the embeddings of node features (i.e., coordinates) and node positional features (i.e., node positions) involve two different aspects of a route solution which are not directly compatible during encoding. They thus proposed the dual-aspect collaborative attention (DACAtt) to learn dual representations for each feature aspect. In this paper, we propose a simple yet powerful Synthesis Attention (Synth-Att) where the attention scores from various types of node feature embeddings can be synthesized to attain a comprehensive representation. It not only has the potential to encode more aspects than DAC-Att, but also reduces the computation costs while reserving competitive performance.

Additionally, we design a diversity enhancement scheme to further ameliorate the performance. The proposed N2S approach was trained through reinforcement learning [Ma et al., 2021b], and we evaluated it on two canonical problems in the PDP family, i.e., the pickup and delivery travelling salesman problem (PDTSP) and its variant with the last-in-firstout constraint (PDTSP-LIFO) to verify our design. Experimental results show that our N2S outperforms the state-ofthe-art, and becomes the first neural method to surpass the well-known LKH3 solver [Helsgaun, 2017] when solving the synthesized PDP instances.

The contributions of the paper are summarized as follows: 1) we present an efficient N2S approach for PDPs, which learns to perform removal and reinsertion of a pickupdelivery node pair automatically; 2) we propose Synth-Att, which allows the vanilla self-attention mechanism to synthesize node relationships from various feature embeddings in a simple way, and achieves superior expressiveness with much fewer computation costs in comparison to DAC-Att; 3) we explore a diversity enhancement scheme, which further makes our N2S the first neural method with almost no domain knowledge to surpass the LKH3 solver on the synthesized PDP instances (when there are no significant distribution shift w.r.t training ones). All these advantages highlight the potential of N2S as a powerful neural solver for new PDP variants without much need for manual trial and error.

2 RelatedWork

Neural Methods for VRPs. We classify recent neural methods into construction and improvement ones. The construction methods, e.g., [Nazari et al., 2018; Joshi et al., 2019; Kool et al., 2018], learn a distribution of selecting nodes to autoregressively build solutions from scratch. Despite being fast, they lack abilities to search (near-)optimal solutions, even if armed with sampling (e.g., [Kool et al., 2018]), local search (e.g., [Kim et al., 2021]), or Monte-Carlo tree search (e.g., [Fu et al., 2021]). Among them, POMO [Kwon et al., 2020] which explored diverse rollouts and data augments is recognized as the best construction method. Differently, improvement methods often hinge on a neighborhood search procedure such as node swap in [Chen and Tian, 2019], ruin-and-repair in [Hottung and Tierney, 2020], and 2-opt in [Wu et al., 2021]. [Ma et al., 2021b] extended the Transformer styled model of [Wu et al., 2021] to Dual-Aspect Collaborative Transformer (DACT), and achieved the stateof-the-art performance, which was also competitive to the hybrid neural methods, e.g., the ones combined with differential evolution [Hottung et al., 2021] and dynamic programming [Kool et al., 2021]. Despite the success of the above methods for CVRP or TSP, they are not verified on the precedence constrained PDPs. Though [Li et al., 2021b] made the first attempt to learn a construction solver for PDPs by introducing the heterogeneous attention to [Kool et al., 2018], the resulting solution qualities are still far from the optimility.

Neighbourhood Search for PDPs. Various heuristics based on neighborhood search have been proposed for PDPs. For PDTSP, [Savelsbergh, 1990] studied the k-interchange neighborhood. A ruin-and-repair neighborhood was later proposed in [Renaud et al., 2000], and further extended in [Renaud et al., 2002] with multiple perturbation methods. Aside from PDTSP, other PDP variants were also studied, normally solved by designing new problem-specific neighborhoods [Parragh et al., 2008], e.g., [Veenstra et al., 2017] proposed five neighborhoods to tackle the PDTSP with handling costs. Among them, the PDTSP-LIFO attracts much attention due to its largely constrained search space. To solve it, [Carrabs et al., 2007] introduced additional neighborhoods such as double-bridge and shake. In [Li et al., 2011], neighborhoods with tree structures were further proposed. Different from the above PDP solvers, the well-known LKH3 solver [Helsgaun, 2017] combines neighborhood restriction strategies to achieve more effective local search, which could solve various VRPs with superior performance. Recently, LKH3 was extended to tackle several PDP variants including PDTSP

(a) δ1 (feasible to both PDPs) (b) δ2 (only feasible to PDTSP)

Figure 2: Two PDP solutions. (a) all goods are on top of the stack at the delivery node; (b) goods from 1+ is blocked by 2+ at 1-.

and PDTSP-LIFO, and delivered comparable performance to [Renaud et al., 2002] and [Li et al., 2011] on PDTSP and PDTSP-LIFO, respectively. Also given the open-sourced nature 2 , we use LKH3 as the benchmark heuristic in this paper.

3 ProblemFormulation

3.1 Notations

We define the studied PDPs over a graph G = (V,E), where nodes in V = P ∪ D ∪ 0 represent locations and edges in E represent trips between locations. With n oneto-one pickup-delivery requests, an PDP instance contains |V | = 2n + 1 different locations, where node 0 is depot, set P = {1+, 2+, ..., n+} contains pickup nodes, and set D = {1-, 2-, ..., n-} contains delivery nodes 3 . Each pickup node i+ has a number of goods to be transported to its delivery node i-. The objective is to find the shortest Hamiltonian cycle to fulfill all requests. In this paper, we consider two representative PDP variants, i.e., PDTSP and PDTSP-LIFO. The solution δ is defined as a cyclic sequence (x0, ..., x2n+1), where x0 and x2n+1 are the depot, and the rest is a permutation of nodes in P ∪ D. For PDTSP, such permutation is under the precedence constraint that requires each pickup i+ to be visited before its delivery i-. For PDTSP-LIFO, the last-in-first-out constraint is further imposed which requires loading and unloading to be executed in the corresponding order. This implies that unloading at a delivery node is allowed if and only if the goods is at the top of the stack. In Figure 2, we present two example solutions with n=3 and |V |=7. The two solutions are both feasible to PDTSP, however, solution δ2 in Figure 2(b) is infeasible to PDTSP-LIFO as the goods from 1+ is NOT at the top of the stack when it needs to be delivered at 1-.

3.2 MDP Formulations

We define the process of solving PDPs by our N2S as a Markov Decision Process M=(S, A, T , R, γ) as follows.

State S. At time t, the state is defined to include: 1) features of the current solution δt, 2) action history, and 3) objective value of the best incumbent solution, i.e.,

where δt is described from two aspects following [Ma et al., 2021b]: l(x) contains 2-dim coordinates of node x (i.e., node features) and pt(x) indicates the index position of x in δt (i.e., node positional feature); H(t,K) stores the most recent K actions at time t if any; and f(·) denotes the objective function and δ∗ t = argminδt′∈{δ0,...,δt} f(δt′ ).

2 http://webhotel4.ruc.dk/∼keld/research/LKH-3/

3 We also refer x0 to the depot, {x1, ..., xn} to the pickup nodes, and {xn+1, ..., x2n} to the delivery nodes from here on.

Action A. With action at ={(i+, i-), (j, k)} where j, k ∈ V \{i+, i-}, the agent removes node pair (i+, i-), and then reinserts node i+ and i- after node j and k, respectively.

State Transition T . We utilize a deterministic transition rule to perform at. An example is illustrated in Figure 1.

Reward R. The reward function is defined as rt =f(δ* t )- min [f (δt+1), f (δ* t )] which is the immediate reduced cost w.r.t. f(δ* t ). The N2S agent aims to maximize the expected total reduced cost w.r.t. δ0 with a discount factor γ < 1.

4 Methodology

4.1 Encoder and Synth-Att

Given state s = {{l(x)}x∈V , {p(x)}x∈V ,H(t,K), f(δ∗ )}, the N2S encoder takes {l(x)}x∈V and {p(x)}x∈V as inputs 4 to learn embeddings for representing the current solution. Following DACT, we first project these raw features into two sets of embeddings, i.e., node feature embeddings (NFEs) {hi}|V | i=0 and positional feature embeddings (PFEs) {gi}|V | i=0 . Different from [Ma et al., 2021b], we treat NFEs as the primary set of embeddings whereas PFEs as auxiliary ones.

NFEs. We define hi as the linear projection of its node features l(xi) for any xi ∈ V with output dimension dh = 128.

PFEs. By extending the absolute positional encoding in the vanilla Transformer [Vaswani et al., 2017], the cyclic positional encoding (CPE) was proposed in [Ma et al., 2021b], which enables Transformer to encode cyclic sequences (as our PDP solutions) more accurately. The PFE gi with output dimension dg = 128 are initialized by CPE as follows,

where superscript d of gi (d) refers to the d-th dimension of gi, and the scalar z(i) as well as the angular frequency ωd are defined in Eq. (3) and Eq. (4), respectively.

According to [Ma et al., 2021b], directly fusing the two sets of embeddings (i.e., hi+gi) may cause undesired noises to the vanilla self-attention. As shown in Figure 3(a), they thus proposed the DAC-Att in ways that each embedding set independently computes attention scores and shares the other aspect to learn dual-aspect representations. Different from it, we propose a simple and generic mechanism by incorporating a multilayer perception (MLP). As shown in Figure 3(b), besides the original self-attention scores (the orange squares), multiple auxiliary attention scores learned from other feature embeddings (the blue squares) are leveraged and fed into an element-wise MLP, which allows it to synthesize heterogeneous attention relationships into comprehensive ones. We call it Synthesis Attention (Synth-Att). It is able to not only leverage more attention scores from various feature embeddings, but also achieve competitive performance to DAC-Att with less computation costs. Below, we present more details.

4H(t,K) is input to decoder and f(δ∗ ) is input to critic network (introduced later). Here we omit time step t for better readability.

Auxiliary Attention Scores. In our N2S, PFEs are used to generate multi-head auxiliary attention scores as follows,

where W Qaux m ∈ Rdg×dq ,W Kaux m ∈ Rdg×dk are trainable matrices for each head m. We set m = 4 and dq = dk = dg/m.

Syn-Att. The Syn-Att is defined as follows,

In specific, it first computes the multi-head self-attention scores αself i,j,m for NFEs based on the trainable matrices W Q m ∈ Rdh×dq andW K m ∈ Rdh×dk for head m using Eq. (7) as

Thereafter, the attention scores αaux and αself are fed into a three-layer MLP with structure (2m× 2m×m) to compute the synthesized multi-head attention scores as follows,

The scores are then normalized to ˜αi,j,m through Softmax, which are further used to calculate the attention values for each head as Eq. (9). Finally, the outputs are given by Eq. (10) with trainable matrix W O ∈ Rmdv×dh (dv = dh/m).

N2S Encoder. We stack L (L=3) encoders, each of which is the same as the Transformer encoder, except that the vanilla multi-head self-attention is replaced with our multihead Synth-Att and we use the same instance normalization layer as [Ma et al., 2021b]. Note that the auxiliary attention scores αaux i,j,m are only computed once and shared among all stacked encoders to reduce the computation costs.

4.2 Decoder

The N2S decoder first adopts the max-pooling layer in [Wu et al., 2021] to aggregate the global representation of all embeddings into each individual one as follows,

Node-Pair Removal Decoder. Given the enhanced embeddings {ˆhi}|V | i=1 and the set H(t,K), the removal decoder outputs a categorical distribution over n requests for removal action. In specific, it first computes a score λi for each xi ∈ V indicating the closeness between node xi and its neighbors as

where pred(xi) and succ(xi) refer to the predecessor and the successor nodes of xi, respectively, and W Q λ ∈ Rdh×dq , W K λ ∈ Rdh×dk . We use multi-head technique to obtain λi,1 to λi,m. Then the decoder aggregates the scores for each pickup-delivery pair (i+, i-) based on a three-layer MLPλ,

where the MLP structure is (2m+4, 32, 32, 1), scalar c(i) counts the frequency of request (i+, i-) being selected for removal in the past K steps, and 1last(t)=i′ is a binary variable indicating whether request i′ was selected at the t-th last step. An activation layer ˆΛ=C · Tanh(˜Λ) is then applied (C=6), followed by Softmax to normalize the distribution which is then used to sample a node pair (i+, i-) as removal action.

Node-Pair Reinsertion Decoder. Given a request (i+, i-) for removal, the reinsertion decoder outputs the joint distribution that reinserts the two nodes back to the solution. We first define two scores µp (xα, xβ) and µs (xα, xβ) for a node xα indicating the degree of preference of accepting a node xβ as its new predecessor and successor nodes, respectively,

where W Qp µ ,W Qs µ ∈ Rdh×dq , and W Kp µ ,W Ks µ ∈ Rdh×dk . Again, we use multiple heads. Based on the scores, the decoder predicts the distribution of reinserting node i+ after node j, and reinserting node i- after node k using MLPµ,

where the MLP structure is (4m, 32, 32, 1). Note that here pred(·) and succ(·) should be considered in the new solution where nodes i+, i- have already been removed. Afterwards, ˆµ = C · Tanh(˜µ) is applied and infeasible choices are masked as -∞ before normalizing by Softmax. Finally, a node pair (j, k), as the reinsertion action, is sampled according to the resulting distribution to indicate the positions of reinserting the node-pair (i-, i+) back to the solution.

Algorithm 1 n-step PPO with CL strategy

Input: policy πθ, critic vϕ, PPO clipping threshold ε, learning rate ηθ, ηϕ, learning rate decay β, epochs E, batches B,

mini-batch κ, training steps Ttrain, CL scalar ρCL 1: for e = 1 to E do 2: for b = 1 to B do 3: Generate training data Db on the fly; 4: Initialize random solutions {δi} to Db; 5: Improve {δi} to {δ′i } via πθ for T = e/ρCL steps; 6: Set initial state s0 based on {δ′i } and Eq. (1); t←0; 7: while t < Ttrain do 8: Get {(st′ ,at′ ,rt′ )}t+n t′=t where at′ ∼πθ(at′ |st′ ); 9: t ← t + n, πold ← πθ, vold ← vϕ; 10: for k = 1 to κ do 11: ˆRt+1 = vϕ(st+1); 12: for t′ ∈ {t, t - 1, ..., t - n} do 13: ˆRt′ ← rt′ + γ ˆRt′+1; 14: ˆ At′ ← ˆRt′ - vϕ(st′ ); 15: end for 16: Compute RL loss JRL(θ) using Eq. (18) and clipped critic loss LBL(ϕ) using Eq. (19); 17: θ ← θ + ηθ∇JRL(θ); 18: ϕ ← ϕ - ηϕ∇LBL(ϕ); 19: end for 20: end while 21: end for 22: ηθ ← βηθ, ηϕ ← βηϕ; 23: end for

4.3 Training Algorithm

As presented in Algorithm 1, the training algorithm for N2S is the same as [Ma et al., 2021b], namely the proximal policy optimization with a curriculum learning strategy. It learns a policy πθ (our N2S) with the help of a critic vϕ as follows.

Critic Network. Given the embeddings {˜h(L) i }|V | i=0 from policy network πθ, the critic network first enhances them by a vanilla multi-head attention layer (with m = 4 heads) to get {yi}|V | i=0 . The enhanced embeddings are then fed into a meanpooling layer in [Wu et al., 2021] to aggregate the global representation of all embeddings into each individual one as

where we use trainable matrices W Local v ,W Global v ∈ Rdh× dh 2 . Lastly, the state value v(st) is output by a four-layer MLP in Eq. (17) with structure (129, 128, 64, 1). Here, f(δ∗ ) is added as an additional input to MLPv.

We train N2S for E epochs and B batches per epoch, where the training dataset Db is randomly generated on the fly with the uniform distribution (line 3). Following [Ma et al., 2021b], we introduce a similar CL strategy (in lines 4 to 6) to further ameliorate the training. In specific, it improves the randomly generated solutions {δi} to {δ′i } by running the current policy πθ for T = e/ρCL steps. The improved solutions

Algorithm 2 N2S-A Inference

Input: Instance I with size |V |, policy πθ, maximum step T 1: for i = 1, ..., ⌊ 12 |V |⌋ do 2: Ii ← I; 3: Ai ←RandomShuffle([flip-x-y, 1-x, 1-y, rotate]); 4: for each augment method j ∈ Ai do 5: ϱj ← RandomConfig(j); 6: Ii ← perform augment j on Ii with config ϱj ; 7: end for 8: end for 9: Solve all instances Ii in parallel with πθ for T steps; 10: return the best solution found among all Ii

{δ′i } with higher quality are then used to initialize the first state s0. In such a way, the hardness level of neighborhood search is gradually increased per epoch following the idea of curriculum learning. The remaining algorithm follows the original design of the n-step PPO, where we present the reinforcement learning loss function in Eq. (18), and the baseline loss function of critic in Eq. (19), respectively. Here, we clip the estimated value ˆ v(st) around the previous value estimates using vclip ϕ (st′) = clip [vϕ(st′), vold(st′) - ε, vold(st′) + ε].

4.4 Diversity Enhancement

To be more resistant to local minima, we further equip our N2S with an augmentation-based inference scheme, which leads to N2S-A in Algorithm 2. Such idea was originally explored in [Kwon et al., 2020] for a neural construction method, and we extend it to an improvement one. The rationale is that an instance I can be transformed into different ones for searching while reserving the same optimal solution, e.g., rotating all locations of nodes by π/2 radian. For an instance of size |V |, our N2S-A performs ⌊ 1 2 |V |⌋ augments, each of which is generated by sequentially applying four preset invariant transformation operations with different orders and different configurations as listed in Table 1. Note that although the mentioned transformations are conducted on instances defined in the Euclidean space, we believe that such an idea has favourable potential to be also exploited in nonEuclidean space, as long as there are proper invariant transformations for coordinates in the target space. Meanwhile, we use K= |V | for training and K′ =⌊ 12 |V |⌋ for inference. This is because we found that when a specific K in H(t,K) is used for training, a smaller K during inference can improve the diversity of the solutions, thus better performance.

5 Evaluation

We design experiments to answer the following questions:

1. How good is the proposed N2S approach against the baselines, including the state-of-the-art neural methods and the strong LKH3 solver? (see Table 3 and Table 4)
2. Can Synth-Att reduce computation costs while achieving competitive performance to DAC-Att? (see Table 5)
3. How crucial are the proposed learnable node-pair removal and node-pair reinsertion decoders for achieving an efficient search? (see Table 6)
4. Can our N2S generalize well to benchmark instances that are different from training ones? (see Table 7)

5.1 Setup

We evaluate N2S on PDTSP and PDTSP-LIFO with three sizes |V | = 21, 51, 101 following the conventions in [Kwon et al., 2020; Ma et al., 2021b], where the node coordinates of instances are randomly and uniformly generated in the unit square [0, 1]×[0, 1]. For our N2S, the initial solution δ0 is sequentially constructed in a random fashion. Our experiments were conducted on a server equipped with 8 RTX 2080 Ti GPU cards and Intel E5-2680 CPU @ 2.4GHz. The training time of our N2S varies with problem sizes, i.e., around 1 day for |V | =21, 3 days for |V | =51, and 7 days for |V | =101, all of which are shorter than the baselines in Table 2.

Hyper-parameters. Our N2S is trained with E = 200 epochs and B = 20 batches per epoch using batch size 600. We set n=5, Ttrain =250 for the n-step PPO with κ=3 minibatch updates and a clip threshold ϵ=0.1. Adam optimizer is used with learning rate ηθ =8×10-5 for πθ, and ηϕ=2×10-5 for vϕ (decayed β=0.985 per epoch). The discount factor γ is set to 0.999 for both PDPs. We clip the gradient norm to be within 0.05, 0.15, 0.35, and set the curriculum learning ρCL to 2, 1.5, 1 for the three problem sizes, respectively.

5.2 Comparison Evaluation

We compare our N2S with the state-of-the-art (SOTA) neural methods and the highly-optimized LKH3 solver.

Regarding the former, we consider the SOTA improvement method DACT [Ma et al., 2021b] 5 and the SOTA construction method Heter-AM [Li et al., 2021b] (specially designed for PDPs). To make a fair comparison with our N2S-A (with the diversity enhancement), we upgrade Heter-AM to HeterPOMO, also given the known superiority of POMO to AM. In specific, we reserve the policy network in Heter-AM as the backbone while adopting the diverse rollouts and the data augmentation techniques in POMO [Kwon et al., 2020] to leverage the advantages of them for the best performance. Each neural baseline is trained using the respective implementation code that is publicly available. For the upgraded Heter-POMO, we adapt and combine the model architecture from the original Heter-AM and the original POMO. The links to their original implementations, approximate training time for the size |V | = 101 , and the used hyper-parameters are presented in Table 2. For other hyper-parameters, we follow the recommendation in their papers. Regarding the latter, LKH3 is a strong heuristic (as reviewed in Section 2) which is widely used as a baseline to benchmark neural methods in recent studies (e.g., [Ma et al., 2021b; Hottung et al., 2021; Li et al., 2021c; Wu et al., 2021]). We report its results with two settings of iterations, i.e., LKH (5k) and LKH (10k).

5 We use the insert decoder which is the best for PDPs.

All baselines are evaluated on a test dataset with 2,000 instances, and we report the metrics of averaged objective values, averaged gaps to LKH3 (10K) and the total solving time. Note that it is hard to perform an absolutely fair time comparison between running Python codes on GPUs (neural methods) and running ANSI C codes on CPUs (LKH solver). Thus we follow the guidelines in [Accorsi et al., 2022] to perform the facilitate comparison that lets each method make full use of the best settings on our machine. In particular, we report the time of LKH3 when running in parallel with 16 CPU cores and the time of each neural method when all 8 GPU cards are available (but do not need to be fully used).

Results on PDTSP. Table 3 shows the results on PDTSP. In the first group, we compare N2S with Heter-AM (greedy and sampling), and DACT. Compared to Heter-AM (5k), our N2S with only 1k steps attains lower gaps with less time for all sizes. Although DACT offers the best gap on PDTSP-21, its performance drops significantly as the problem size increases, partly because its decoder is less efficient than our node-pair removal and reinsertion ones when tackling largerscale problems. Instead, our N2S achieves higher performance and consistently dominates DACT in terms of both the gaps and the time on PDTSP-51 and PDTSP-101. In the second group, our augmented N2S-A is compared to the upgraded Heter-POMOmethod with three variants 9 . It is shown that even with only 1k steps, our N2S-A attains a significantly smaller gap than all three Heter-POMO variants, by almost an order of magnitude. Although Heter-POMO-A (gr.) tends to be competitive with fast speed, the gap is hard to be further reduced by increasing inference time if we refer to HeterPOMO-A (3k). Moreover, our N2S-A (2k) keeps abreast of, or even slightly exceeds the strong LKH3 solver, achieving gaps of -0.03% and -0.01% with less time on PDTSP-51 and PDTSP-101, respectively. Those gaps are further reduced to -0.04% and -0.20% with more steps, i.e., 3k.

6 https://github.com/Demon0312/Heterogeneous-AttentionsPDP-DRL

7 https://github.com/yd-kwon/POMO

8 https://github.com/yining043/VRP-DACT

9 Heter-POMO (gr.), Heter-POMO-A (gr.), Heter-POMO-A (3k) refer to: greedily generate |V | solution; greedily generate |V | solutions with 8 augments; sample |V |×3k solutions with 8 augments.

Results on PDTSP-LIFO. In Table 4, we report the results on PDTSP-LIFO. Due to the more constrained search space, DACT failed to work well (see Table 6). Therefore, we mainly compare our N2S-A with Heter-POMO (the best neural baseline in Table 3) and the LKH3 solver. As exhibited, the advantages of neural methods over LKH3 have been further enhanced on this harder problem, where our approach consistently outperforms Heter-POMO for all sizes. Compared to LKH3, our N2S-A with 3k steps presents superior performance again, and attains gaps of -0.64% and -1.47% on PDTSP-LIFO-51 and PDTSP-LIFO-101, respectively.

5.3 Ablation Evaluation

Effects of Different Encoding Methods. We replace the proposed Synth-Att in our N2S encoder with vanilla-Att and DAC-Att, respectively. Using only one GPU card, we report the number of model parameters, time, and gaps for solving 2,000 PDTSP-51 instances with 3k steps in Table 5. As exhibited, Synth-Att attains slightly smaller gaps than DAC-Att with much fewer computation costs, and considerably lower gaps than vanilla-Att with slight extra computation costs.

Effects of Different Decoding Methods. To highlight the desirability of our two decoders, we replace the trainable decoders with hand-crafted ones (i.e., random and the ϵ-greedy with ϵ = 0.1) while ensuring feasibility. We gather the results on PDTSP-51 and PDTSP-LIFO-51 with 3k steps in Table 6 where mark '✓' means our proposed decoder is used (the one without augments) and mark '✗' means the decoder in the parentheses is used instead. As revealed, the methods of retaining at least one trainable decoder always attain much lower gaps than the ones equipped with only hand-crafted decoders. The method with both trainable decoders (i.e., N2S)

achieves the best performance. We also notice that DACT fail to solve PDTSP-LIFO well. This might be because its decoder considers removing and reinserting only one node instead of a pickup-delivery node pair in each action, which is a serious limitation for more constrained PDPs. For example, on PDTSP-LIFO, over 92% of the action space need to be masked for DACT, which leads to extremely low efficiency.

5.4 Generalization Evaluation

We further evaluate our N2S on benchmark instances, including all the ones from [Renaud et al., 2002] for PDTSP and the ones with size |V | ≤201 from [Carrabs et al., 2007] for PDTSP-LIFO, which are largely different from our training ones, e.g., different sizes (i.e., 200 nodes) as shown in Figure 4 and different node distributions as shown in Figure 5. In Table 7, we report the best and the average gaps (with 10 independent runs) achieved by N2S-A and neural baseline Heter-POMO-A w.r.t. optimal solutions for PDTSP, or heuristic baseline B1 [Carrabs et al., 2007] and B2 [Li et al., 2011] for PDTSP-LIFO. It can be seen that our N2S significantly outstrips Heter-POMO in all cases. Without re-training or leveraging other techniques, this is fairly hard to achieve because machine learning often suffers from a mediocre outof-distribution zero-shot generalization [Zhou et al., 2021; Li et al., 2021a]. The results also imply slight inferiority of our N2S to the LKH3 solver, given that LHK3 reports similar performance to the B2 baseline on those instances [Helsgaun, 2017]. Accordingly, we will focus on improving the out-ofdistribution generalization for our N2S in future.

6 Conclusion

We present an efficient N2S approach for PDPs. It utilizes a novel Synth-Att to synthesize node relationships from various types of solution features and exploits the node-pair removal and reinsertion decoders to tackle the precedence constraint. Extensive experiments on PDTSP and PDTSP-LIFO verified our design, where N2S achieves state-of-the-art performance among existing neural methods. Further equipped with a diversity enhancement scheme, it even becomes the first neural method to surpass LKH3 on synthesized PDP instances. In future, we will 1) deploy N2S as a low-level agent in the hierarchical framework of [Ma et al., 2021a] for dynamic PDPs, 2) combine N2S with a similar divide-and-conquer strategy in [Li et al., 2021c] for much larger-scale instances, and 3) enhance the out-of-distribution generalization to exceed LKH3 on instances of arbitrary distributions.

Acknowledgments

This work was supported by the IEO Decentralised GAP project of Continuous Last-Mile Logistics (CLML) at SIMTech under Grant I22D1AG003. It was also supported in part by the National Natural Science Foundation of China under Grant 62102228, and in part by Guangdong Regional Joint Fund for Basic and Applied Research under Grant 2021B1515120078.