A case study of the profit-maximizing multi-vehicle pickup and delivery selection problem for the road networks with the integratable nodes

This paper is a study of an application-based model in profit-maximizing multi-vehicle pickup and delivery selection problem (PPDSP). The graph-theoretic model proposed by existing studies of PPDSP is based on transport requests to define the corresponding nodes (i.e., each request corresponds to a pickup node and a delivery node). In practice, however, there are probably multiple requests coming from or going to an identical location. Considering the road networks with the integratable nodes as above, we define a new model based on the integrated nodes for the corresponding PPDSP and propose a novel mixed-integer formulation. In comparative experiments with the existing formulation, as the number of integratable nodes increases, our method has a clear advantage in terms of the number of variables as well as the number of constraints required in the generated instances, and the accuracy of the optimized solution obtained within a given time.


Introduction
Recently, the profit-maximizing multi-vehicle pickup and delivery selection problem (PPDSP) has gained a lot of attention in the field of practical transportation and logistics Liu et al. (2018); Riedler & Raidl (2018); Asghari & Mirzapour Al-e-hashem (2020); Huang et al. (2020). This problem was first proposed in Qiu et al. (2017), which involves three classical problem models: routing optimization, pickup and delivery, and selective pickup (a.k.a. knapsack ). Solving this problem quickly and optimally can both help improve the operational efficiency of the carriers and contributes to more eco-friendly transportation.
To the best of our knowledge, in the graph-theoretic models constructed in the existing studies Ting et al. (2017); Gansterer et al. (2017); Al-Chami et al. (2018); Ahmadi-Javid et al.   Figure 1c can be regarded as a feasible solution of the request-based model, but cannot be corresponded to in the location-based model.
(2018), the definitions of nodes are based on the pickup and delivery location from the requests. This means that one request needs to correspond to two nodes. However, in application scenarios, there are often plural requests coming from or to arrive at the same location. The number of variables to be required in the existing mixed-integer formulation heavily depends on the number of nodes in the model, and since PPDSP is an N P-hard problem, its computational complexity is exponential with respect to the number of variables.
The motivation of this study is to provide a more reasonable and effective mathematical model for practical logistics and transportation problems. Considering the road networks with the integratable nodes as above, we proposes a new concise model in which the definition of nodes is based on the locations rather than the requests, and give a novel mixed-integer programming formulation that reduces the number of required variables. It is necessary to claim that the solution set of our proposed location-based method is a subset of the solution set of the request-based method. This is because the integrated node can only be passed at most once under the Hamiltonian cycle constraint, where "at most" is due to the objective function of profit-maximizing (i.e., it is possible that any node will not be traversed). In Figure 1, we take the delivery of single truck as an example, and assume that the truck is transported according to the route given in Figure 1a as a feasible solution in this scenario, where the pickup coordinates of request 1 (abbr. Req. 1 pickup coord.) and Req. 2 pickup coord. are the identical locations. We can integrate them as location 1 as shown in Figure 1b, and after integrating all the same coordinates into one location respectively, our model becomes more concise. This integration changes the nodes of the model from repeatable coordinates based on the request to the unique locations, which results in each location can only be visited once. Therefore, the location-based model can also correspond to the feasible solutions in Figure 1a, but not to some of the feasible solutions corresponded to the request-based model, such as the solution given in Figure 1c. Whereas in real logistics application scenarios, making multiple retraces to the same location is rare (e.g. insufficient vehicle capacity, goods cannot be mixed, etc.), and also non-efficient. For such locations with a large number of requests or high loading volume, it is a more common strategy to assign multiple vehicles to them. Therefore, we argue that the location-based model can improve the optimization efficiency although it reduces the range of feasible solutions.

Preliminaries
Given a directed graph G = (V, E), where V = {0, 1, 2, . . . , |V |} is the set of nodes representing the location points (0 is depot) and E is the set of arc, denoted as > od, and given a set of trucks T = {1, 2, . . . , |T |}, we define a type of Boolean variables x t od that is equal to 1 if the truck t passes through the arc > od and 0 otherwise. We denote the load capacity of truck t as c t , and the cost of the truck t traversing the arc > od as l t od , where t ∈ T and o, d ∈ V .
Let R = {1, 2, . . . , |R|} be a set of requests. Each request r (r ∈ R) is considered as a tuple r = w r , q r , f (r), g(r) , where • w r is the the payment that can be received for completing the shipping of request r; • q r is the volume of request r; • f (r) is the loading point of request r; • g(r) is the unloading point of request r, and f : R → V \ {0} (resp. g : R → V \ {0}) is a function for mapping the loading (resp. unloading) point of requests r. We also define another type of Boolean variables y t r that is equal to 1 if request r is allocated to truck t and 0 otherwise, where r ∈ R and t ∈ T . Besides, we denote the number of location points visited by truck t before it reaches Definition 1 (Delivery of Truck t). Delivery of truck t is denoted by D t = r∈R {r | y t r = 1}, where t ∈ T .
Definition 2 (Route of Truck t). Route of truck t is denoted by } as the set of location points visited and departed exactly once by truck t, where the predecessor and successor are the same node is not counted (e.g., even if > vv | x t vv = 1, neither this time can be included in the number of visits or departures of truck t to/from location v.); • Ensure that no subtour exists in S t . 1 Loading Before Unloading ∀i ∈ D t , u t f (i) < u t g(i) .
Capacity Limitation At any time, the total volume of cargo carried by truck t cannot exceed its capacity c t .
Definition 3 (Delivery Routing Solution). Delivery routing solution DS = t∈T {(D t , S t )} that is a partition of R into disjoint and contained D t with the corresponding S t : Denote the set of all possible delivery routing solutions as Π(R, T ).
Definition 4 (Profit-Cost Function). A profit-cost function assigns a real-valued profit to every D t : w : D t → R and a cost to every S t : l : S t → R, where w(D t ) = r∈Dt w r and l(S t ) = > od∈St l t od . For any delivery routing solution DS ∈ Π(R, T ), the value of DS is calculated by In general, a delivery routing solution DS that considers only maximizing profits or minimizing costs is not necessarily an optimal DS. Therefore, we have to find the optimal delivery routing solution that maximizes the sum of the values of profit-cost functions. We define a delivery routing problem in profit-cost function as follows.
Definition 5 (PPDSP). For a set of requests and trucks (R, T ), a profit-maximizing multivehicle pickup and delivery selection problem (PPDSP) is to find the optimal delivery routing solution DS * such that Here we show an example of PPDSP.
Example 1. Assume that two trucks T = {t 1 , t 2 } are responsible for three requests R = {r 1 , r 2 , r 3 } with the following conditions.
• The information about the requests: • The information about the trucks: 4 -The capacities of the trucks c t 1 = 6 and c t 2 = 3.
-The cost matrices of each truck through each arc (δ is depot).
In this example, the optimal delivery routing solution DS * for the given PPDSP is cδ})}, and its value is 11.

Problem Formulation
In this section, we present the following mixed-integer programming (MIP) formulation of the location-based model for PPDSP and prove its correctness as well as the space complexity 5 of the generated problem.

Correctness
The objective function in Eq. (1) maximizes the profit-cost function for all possible delivery routing solutions. Eq. (2) introduces two types of Boolean variables x t od and y t r . Eq. (3) guarantees that, each request can be assigned to at most one truck. Proof. According to Eq. (4) (resp. Eq. (5)), if request r is assigned to truck t, then truck t reaches the pickup (resp. dropoff) point of r at least once via a point that is not the pickup (resp. dropoff) point of r. Eq. (6) ensures that the number of visits and the number of departures of a truck at any location must be equal. Eq. (7) restricts any truck departing from a location to at most one other location. It is clear from the above that, all nodes of r∈Dt {f (r), g(r)} are ensured to be visited and departed by truck t exactly once.
Eq. (8), a canonical MTZ subtour elimination constraint Miller et al. (1960), resrticts the integer variables u t v , where v ∈ V \ {0}, such that they form an ascending series to represent the order that truck t arrives at each location v. Eq. (12) gives the domain of such variables u t v . Furthermore, Eq. (8) associated with Eqs. (4)-(7) also enforces that depot 0 must be visited and departed by truck t exactly once if D t = ∅.
Eq. (9) ensures that, for any request r, the arrival of its pickup point must precede the arrival of its dropoff point by truck t if y t r = 1. Theorem 2. For PPDSP, the capacity constraint can be guaranteed by the Eqs. (10) and (11). (10), whose domains are between 0 and c t (i.e., the capacity of truck t) as given in Eq. (11). Such a variable h t v can be considered as the loading volume of truck t when it departs at location v. We define the amount of change in the loading of truck t at location v as Γ, which equals the loaded amount at v (i.e., r∈f −1 (v) (q r · y t r )) minus the unloaded amount at v (i.e., r∈g −1 (v) (q r · y t r )). The main part of Eq. (10) guarantees that h t d − h t o is exactly equal to Γ when x t od = 1, which is written as a big-M linear inequality, where we specify M as c t + r∈R q r .

Space Complexity
Consider that the number of trucks is m (i.e., |T | = m) and the number of requests is n (i.e., |R| = n). Assume that the average repetition rate of the same locations is k, where k ≥ 1, then we can estimate the number of unique locations as 2nk −1 , which is also the number of integrated nodes (viz., |V \ {0}| = 2nk −1 ). Proof. The number of Boolean variables x t od (resp. y t r ) required in proposed formulation is m(1 + 2nk −1 ) 2 (resp. mn); while both the number of required integer variables u t d and h t d are m(1 + 2nk −1 ). In Table 1, we list the bounded numbers of linear (in)equations that correspond to Eqs. (3)-(10) involved in the proposed PPDSP formulation. Therefore, the number of required variable as well as the number of corresponded linear (in)equations are of Θ(mn 2 k −2 ).

Experiment
We compare the performance of a MIP optimizer in solving randomly generated PPDSP instances based on the proposed formulation (i.e., location-based model), and based on the existing formulation (i.e., request-based model), respectively. Please refer to Appendix A for the specific description of the existing formulation used for the comparative experiments. 7

Instances Generation
The directed graph informations for generating instances are set based on the samples of TSPLIB benchmark with displaying data as the coordinates of nodes. 2 We denote the number of nodes contained in the selected sample as |V |, and set the first of these nodes to be the depot node. The pickup and dropoff points of each request are chosen from the non-depot nodes. Therefore, given an average repetition rate k for the same locations, we need to repeatably select n pairs of non-depot nodes (i.e., a total of 2n repeatable non-depot nodes) to correspond to n requests, where n = round( k(|V |−1) 2 ). In Algorithm 1, we construct the repeaTimeList of length |V | − 1 to mark the number of times of each non-depot node being selected. In order to ensure that each non-depot node is selected at least once, each value on such list is initialized to 1. We continuously generate a random index of repeaTimeList and add one to the value corresponding to the generated index (i.e., the number of times of the non-depot node being selected increases by one) until the sum of these numbers reaches 2n. The pickup and dropoff nodes for each request are randomly paired up in Algorithm 2. We first rewrite repeaTimeList as a list of length 2n, shuffList, which consists of all selected non-depot nodes, where the number of repetitions indicates the number of times they are selected. For example, we have shuffList = [0, 0, 0, 1, 2, 2] for repeaTimeList = [3, 1, 2]. Next, we shuffle shuffList and clear pairList, which is used to store pairs of nodes indicating the pickup and dropoff points of all requests. We then divide shuffList into n pairs in order.
2 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/ 8 If two nodes of any pair are identical (i.e., the pickup and dropoff are the same point), or the pair with considering order is already contained in pairList, then we are back to Line 4. Otherwise, we append the eligible pair of nodes into pairList until the last pair is added.
Algorithm 2: Randomly pair up the pickup and dropoff nodes for each request if i = n − 1 then reshuffle ← false 12 return pairList Since we intend to generate instances corresponding to different k for each selected sample of TSPLIB benchmark, and the number of requests n gets smaller as k decreases, we need to produce a decrementable list of node pairs such that we can remove some of the node pairs to correspond to smaller k, but the remaining part of the node pairs still contains all locations in the sample other than that as depot (i.e., they still need to be selected at least once). Therefore, in Algorithm 3, we assume that a sorted list of node pairs pops the end elements out of it one by one, yet it is always guaranteed that every non-depot node is selected at least once. We insert the node pair in which both nodes are currently selected only once into the frontmost of head in Line 7, and then append the node pair in which one of them is selected only once into head in Line 11. Then we keep inserting the node pair in which the sum of the repetitions of the two nodes is currently maximum into the frontmost of tail in Line 22. Note that such insertions and appends require popping the corresponding node pair out of pairList and updating L. We merge head and tail to obtain a sorted list of node pairs sortedPairs at the end of Algorithm 3. Finally, we can generate a list containing n requests for each selected sample of TSPLIB benchmark, as shown in Algorithm 4.
In addition, we generate data for the three types of trucks recursively in the order of maximum load capacity of 25, 20, and 15 until the number of generated trucks reaches m. And these three types of trucks correspond to their respective cost coefficients for each traversing arc of 1.2, 1, and 0.8. For example, for a truck t of load capacity 15, l t od , i.e., the Algorithm 3: Sort the node pairs by the sum of the repeat times of each node in the node pair Input: repeaTimeList, pairList We generate a total of 3 × 5 × 5 = 75 PPDSP instances for our comparative experiments based on the proposed formulation and on the existing formulation, respectively, according to the following parameter settings.
• The selected TSPLIB samples are {burma14, ulysses16, ulysses22 }; 3 Distance( > od) is the Euclidean distance between the coordinates of node o and the coordinates of node d.

Experimental Settings
All experiments are performed on an Apple M1 Pro chip, using the Ubuntu 18.04.6 LTS operating system via the Podman virtual machine with 19 GB of allocated memory. A single CPU core is used for each experiment. We implemented the instance generators for both the proposed formulation and the existing formulation by using Python 3. Each generated problem instance is solved by the MIP optimizer-Cplex of version 20.1.0.0 CPLEX (2020) within 3,600 CPU seconds time limit. Source code for our experiments is available at https://github.com/ReprodSuplem/PPDSP.

Results
Tables 2-4 show the performance of the existing formulation-based method and our proposed formulation-based method in terms of the number of generated variables (#Var.), the number of generated constraints (#Con.), and the optimal values (Opt.) for the various average repetition rates of the same locations (k) and the different numbers of trucks (m), corresponding to TSPLIB samples burma14, ulysses16 and ulysses22, respectively. Each cell recording the left and right values corresponds to a comparison item, where the left value refers to the performance of the item based on the existing formulation-based method, while the right value corresponds to the performance of the item based on the proposed formulation-based method. We compare the performance of these two methods in such cells and put the values of the dominant side in bold.
We can see that either the number of variables or the number of constraints generated by our proposed method is proportional to m, while neither the number of variables nor the number of constraints generated by our proposed method increases significantly as k becomes larger. Such a result is consistent with Theorem 3. In contrast, although the number of variables and the number of constraints produced by the existing method is Table 2: Comparison of the existing formulation-based method with the proposed formulation-based method for TSPLIB sample burma14 in terms of the number of generated variables, the number of generated constraints, and the optimal values. m item k = 1 k = 1.5 k = 2 k = 2.5 k = 3 (n = 7) (n = 10) (n = 13) (n = 16) (n = 20) also proportional to m, with k becoming larger, exponential increases in both of them are observed. Furthermore, it is interesting to note that even when k = 1, the proposed method generates fewer variables than the existing method. As k increases, the optimal value obtained by the proposed method increases within the time limit; while there is no such trend in the optimal values obtained by the existing method. In addition, as m grows up, in theory, the upper bound of the optimal value cannot be smaller. However, the optimal values obtained by the proposed method and the existing method do not maintain the trend of increasing, instead, they both have inflection points. This is due to the fact that the search space of the problem becomes huge, which makes it inefficient for the solver to update the optimized solution. There are even cases on the existing method where the initial solution is not obtained until the end of time.
Last but not least, we can clearly see that for those problems with larger k (i.e., more nodes that can be integrated), the method based on our proposed formulation generates fewer variables and constraints, as well as achieves larger optimal values, than the method based on the existing formulation.

Conclusion
In this paper, we revisit PPDSP on road networks with the integratable nodes. For such application scenarios, we define a location-based graph-theoretic model and give the corresponding MIP formulation. We prove the correctness of this formulation as well as analyze its space complexity. We compare the proposed method with the existing formulation-based method. The experimental results show that, for the instances with more integratable nodes, our method has a significant advantage over the existing method in terms of the generated problem size, and the optimized values.

A. Existing Formulation
Here we provide a detailed description of the existing formulation proposed in Qiu et al. (2017) that were used to perform the comparative experiments in Section 4.
As we stated before, in order to compare the most essential differences between the request-based model and the location-based model, the constraints on time windows are not considered in this study. However, because in the existing formulation, the time windows constraints also act as the subtour elimination constraint. For this reason, we replace those constraints on time windows in the existing formulation with the same MTZ subtour elimination constraint as in our proposed formulation to ensure the fairness of the comparative experiments.
In the existing formulation, n requests are represented as a directed graph G = (V, E), whose node set V = {0} ∪ P ∪ D ∪ {2n + 1}, where P = {1, 2, . . . , n} is the set of pickup nodes, D = {n + 1, n + 2, . . . , 2n} is the set of dropoff nodes, and node 0 and 2n + 1 indicate two depots. 4 A request r is formalized by using the combination (r, r + n) to refer to the pickup and dropoff nodes, and each node v (v ∈ P ∪ D) corresponds to a change amount q v of the truck's load, where q r and q r+n respectively correspond to the loading and unloading of request r, q r > 0, q r+n < 0, q r = |q r+n | and r ∈ {1, 2, . . . , n}. In addition, we particularly specify q 0 = q 2n+1 = 0 for the two depot nodes. The definitions and notations of the other variables are consistent with that previously described in this paper. The existing formulation used for comparative experiments are as follows.
It should be noted that we found three incorrectnesses in the existing formulations in Qiu et al. (2017), which correspond to Eqs. (A.3), (A.4) and (A.12) in the above formulation, and we have rectified each of them.