Multi-start heuristic approaches for one-to-one pickup and delivery problems with shortest-path transport along real-life paths

The One-to-one Pickup and Delivery Problem with Shortest-path Transport along Real-life Paths (OPDPSTRP) is presented in this paper. It is a variation of the One-to-one Pickup and Delivery Problem (OPDP), which is common in daily life, such as the Passenger Train Operation Plans (PTOP) and partial Taxi-sharing Problem. Unlike the classical OPDP, in the OPDPSTRP, (1) each demand must be transported along the shortest path according to passengers/shippers requirements, and (2) each vehicle should travel along a real-life path. First, six route structure rules are proposed for the OPDPSTRP, and a kind of Mixed-Integer Programming (MIP) models is formulated for it. Second, A Variable Neighborhood Descent (VND), a Variable Neighborhood Research (VNS), a Multi-Start VND (MS_VND) and a Multi-Start VNS (MS_VNS) with five neighborhood operators has been developed to solve the problem. Finally, The Gurobi solver, the VND, the VNS, the MS_VND and the MS_VNS have been compared with each other by 84 random instances partitioned in small size connected graphs, medium size connected graphs and large size connected graphs. From the test results we found that solutions generated by these approaches are often comparable with those found by the Gurobi solver, and the solutions found by these approaches are better than the solutions found by the Gurobi solver when solving instances with larger numbers of demands. In almost all instances, the MS_VND significantly outperforms the VND and the VNS in terms of solution quality, and outperforms the MS_VNS both in terms of solution quality and CPU time. In the instances with large numbers of demands, the MS_VND is still able to generate good feasible solutions in a reasonable CPU time, which is of vital practical significance for real-life instances.


OPDP, OPDPST and OPDPSTRP
Most classical OPDPs are studied in complete graphs, and pickup points must be visited prior to delivery points (e.g. Fleischmann et al. [21], Cornuéjols et al. [22], Aragão et al. [23], Letchford et al. [24], Li et al. [25], Mahmoudi et al. [26]). The classical OPDP can be easily formulated as a Mixed-Integer Program (MIP), such as those reported in Pérez et al. [17], Sahin et al. [18] and Soysal et al. [19]. In a forthcoming article, we studied a new kind of OPDP, named One-to-one Pickup and Delivery Problems with Shortest-path Transport (OPDPST), in which each pd-pair must be transported along the shortest path, and a new kind of modeling method was proposed for the OPDPST according its new route constructions. The OPDPSTRP is a kind of the OPDPST, in which vehicles should travel along real-life paths in connected graphs, such as Train-operation Plan in China and partial Taxi-sharing Problem and so on. Unlike the classical OPDP in a complete graph, the OPDPST and the OPDPSTRP studied in this paper is described by connected graphs, since numbers of vehicle stops need to be considered.
Vehicles should travel along paths both in the classical OPDP and the OPDPST and the OPDPSTRP, some research conducted on the route structure of the classical OPDP can provide some reference for the OPDPSTRP, although travelling along paths in a complete graph doesn't means travelling along paths in the corresponding real-life connected graph. Furuhata et al. [27] listed four kinds of ride-sharing patterns for the Ride-sharing Problem. Letchford et al. [24] pointed out that the cheapest path is not always the quickest path, and a comparison of multiple paths between each two points was necessary. Muelas et al. [28] proposed a method for relocating a pd-pair by considering four cases, and the shortest one was chosen as the optimal routing scheme in each local search move. Lin et al. [29] and Wang et al. [30] studied the Ride-sharing Problem (a kind of OPDP) in real-life networks. However, it is not a necessary requirement to transport pd-pairs through the shortest route in these papers.
Additionally, as in the OPDPSTRP, each vehicle starts at its location (regarded as a depot) and ends at the final delivery point of the contents transported by the vehicle, so it can be considered as a multi-depot (vehicles) problem. Most OPDP research is based on single depot, such as that reviewed by Psaraftis [12], [13] (1983), Desrosiers et al. [31], Lin et al. [29], Li et al. [32], Letchford et al. [24], Urra et al. [33], Muelas et al. [28], Qiu et al. [34] and Li et al. [25]. There is also some OPDP research based on multiple depots (vehicles), which is mainly concerned with the Taxi-sharing Problem and Ride-sharing Problem. For example, there is a starting point and an ending point for each vehicle in Hosni et al. [35], Ma et al. [36], Sawaya et al. [37] and Santos et al. [9], whilst only the starting point is considered for each vehicle in Detti et al. [38].
Generally, there is far more research on classical OPDP, and little research focuses directly on the OPDPSTRP proposed in this paper. Cordeau et al. [39] presented eight kinds of local search moves for OPDP: Couple-exchange, Blockexchange, Relocate-couple, Relocate-block, Multi-relocate, 2-opt-L, Double-bridge and Shake. Ribeiro et al. [40] and Ropke et al. [41] modified three large neighborhood removal heuristics and two large neighborhood insertion heuristics from Shaw P. [42], [43] and Potvin et al. [44] for OPDP. Additionally, studies of Grimault et al. [45] and Ho et al. [7] show that solution feasibility of the OPDP is an important issue to the neighborhoods efficiency for the algorithm.

Neighborhood and algorithm for OPDP
Rodríguez-Martín et al. [46] solved the OPDPTSP by the Greedy Randomized Adaptive Search Procedures (GRASP) and the VND. Sahin et al. [18] proposed an efficient heuristic that combines the strengths of Tabu Search and Simulated Annealing for the OPDPSD. A Iterated Local Search (ILS) is proposed by Li et al. [47]. Ho et al. [7] classified the solution methods for the DARP (an important category of the OPDP). Factorovich et al. [48] propose a combination of cutting planes to find feasible solutions for the pickup and delivery problem with incompatibility constraints. Additionally, some Local Search (LS) meta-heuristics studied for the PDP can also be used for reference. Adaptive Large Neighborhood Search (ALNS) proposed for the PDP by Ropke et al. [41] and Li et al. [25]. Berbeglia et al. [8], [9] reviewed the algorithms for the static and the dynamic PDP. There are also some exact methods proposed for the OPDP, Pérez et al. [17] solve two mixed integer linear programming models of the OPDPTSP by the Cplex solver. Soysal et al. [19] proposed a mixed integer programming model for the green OPDP, and solved it by the Cplex solver.
Generally, research for the OPDP is mainly based on the heuristic techniques, such as the VNS, the VND, the ILS, the ALNS, the hybrid algorithms and so on, because exact methods are often not the most effective way to solve large OPDP.

Problem definition
In order to define the proposed the OPDPSTRP in mathematical terms, we specify an connected graph G = (N,E,P,K), in which N = {1,. . .,n} for vertexes, E = {1,. . .,e} for edges, P = {1,. . .,p} for pd-pairs, and K = {1,. . .,m} for vehicles. Each pd-pair i with demand q i yields revenue π i . Each vehicle k2K has a maximum capacity Q k and a fixed cost vc k . The transportation cost per unit length of vehicle k is tc k . Each vehicle k has a stop cost sc k n at node n. The system also obeys the following assumptions.
1. Each pd-pair attribute is different (Perez et al. [17], Psaraftis et al. [49], Miao et al. [50]) and cannot be split (unlike Mitra et al. [51] and Nowak [52], Xia et al. [53]), and must also be transported through the shortest path according to passengers requirements, with pickup point being visited prior to delivery point.
2. Each vehicle must travel along a real-life path beginning with the first pickup point and end at the last delivery point, namely each point cannot be accessed multiple times by a vehicle, a common practice in the Passenger Train Operation Plans and so on. The model with this constraint always can obtain close solutions to the model without this constraint but with less computing time than it, which is elaborated in another initial manuscript, the major results can be found in S3 Appendix.

Constants and variables
The constants and variables used in this paper are listed in Table 1.

Mathematical model
The model of the OPDPSTRP is formulated in this section, and the route structure of the OPDPSTRP will be studied for it in Section 3.4.
y k e 2 f0; 1g sn k n 2 f0; 1g The objective function (1) maximizes the total profit, in which ;j is the total fixed cost of using vehicles, X k2K X e2E tc k � le e �y k e is the total travel cost due to carry lengths, tc k � lc i;j � x k i;j is the total traveling cost due to connection lengths and additional lengths, and X k2K X n2N sc k n � sn k n is the total cost of stop nodes. Constraints (2), (3), (4), (5) and (6) determine the order between pd-pairs/vehicle i and j. Constraints (7) ensure that each pd-pair i is transported no more than once. Constraint (8) ensures that each vehicle is not over-loading. Constraint (9) determines whether the vehicle stops at node n or not. Constraint (10) ensures that the number of stops for each route (not including the depot/vehicle) does not exceed M. Constraint (11) determines whether edge e is traveled along or not. Constraints (12) ensure that each route with pd-pairs is assigned to one vehicle. Constraint (13) ensures that the length of each route (not including the depot/vehicle) is not longer than D. Constraints (14), (15), (16) and (17) introduce the decision variables.

Route structure of the OPDPSTRP
The route structure of the OPDPSTRP will be studied for the model. 3.4.1 Feasibility of route structure and connection relationship between two pd-pairs or vehicle/pd-pair. Definition 1: In a real-life connected graph, if all pd-pairs are transported through the shortest path in a route i starting with the first pickup point, and the route is a path, then it is defined that route i is Route-Structure-Feasible (RSF for short).
Definition 2: If a RSF route stems from inserting pd-pair i into route j, then it is defined that pd-pair i can be inserted into route j according route structure. pd_R_rs_judge i,j is defined as the route structure feasibility judgement parameter of inserting (RSFJPI for short, 0: feasible; 1: unfeasible.) pd-pair i into route j. [pd_R_rs_judge i,j ] is defined as the route structure feasibility judgement matrix of inserting (RSFJMI for short).
Definition 3: If a RSF route is a combination of pd-pair i and pd-pair j, then it is defined that pd-pair i can be combined with pd-pair j according route structure. pd_combine_rs_judge i,j is defined as route structure feasibility judgement parameter of combining (RSFJPC for short, 0: feasible; 1: unfeasible) pd-pair i with j. [pd_combine_rs_judge i,j ] is defined as the route structure feasibility judgement matrix of combining (RSFJMC for short). It is assumed that each pd-pair can combine with any vehicle, that is pd_combine_rs_judge p+1,j = 1(8j = 1,. . .,p).
Definition 4: In a RSF route, if pd-pair j can be picked up not prior to vehicle/pd-pair i, then it is defined that pd-pair j can connect to vehicle/pd-pair i. The parameter connect_to_ judge i,j (8i = 1,. . .,p+1;j = 1,. . .p) is defined as the judgment parameter for pd-pair j connecting to vehicle/pd-pair i. It is assumed that each pd-pair can connect to any vehicle, that is con-nect_to_judge p+1,j = 1(8j = 1,. . .,p). One vehicle cannot connect to another vehicle.
Definition 5: If pd-pair i can connect to vehicle/pd-pair i (namely connect_to_judge i,j = 1), and pd-pair j can be delivered not prior to vehicle/pd-pair i, then it is defined that pd-pair j can connect after vehicle/pd-pair i. The parameter connect_after_judge i,j (8i = 1,. . .,p+1;j = 1,. . .p) is defined as the judgment parameter for pd-pair j connecting after vehicle/pd-pair i. It is assumed that each pd-pair can connect after any vehicle, that is, connect_after_judge p+1,j =1 (8i = 1,. . .,p).
Definition 6: The nodes traveled by vehicle k in the route combined by pd-pair i and pdpair j can be classified into stop nodes (stopped at by vehicle k) and pass nodes (passed but not stopped at by vehicle k). p i and d i are pickup point and delivery point of pd-pair i correspondingly. sn_od i,n is defined as the judgment parameter of whether pd-pair i picks-up/delivers at node n.
Definition 7: The sections in the route combined with pd-pair i and pd-pair j include weighting sections (sections traveled by a vehicle with pd-pairs), and connecting sections (sections traveled by a vehicle without pd-pairs). The constants le e and lc i,j are defined as their lengths. Table 2 shows the lengths of sections in the paths (beginning with the first pickup point) constructed by pd-pair i and pd-pair j in Fig 2. 3.4.2 Rules of route construction with more than two pd-pairs/vehicles. Since x k i;j is defined as the decision variable of whether pd-pair j connects to pd-pair i or not, and connec-t_after_judge i,j = 1 means pd-pair j can connect after pd-pair i, so x k i;j � connect a fter j udge i;j ¼ 1 means pd-pair j connects after pd-pair i successfully (Definition 4 and Definition 5). According to the above research, in a route combined with more than two pd-pairs/vehicle traveled by vehicle k, the lengths of the weighting sections and connecting sections can be formulated as X e2E le e � y k e and X i2P X j2P lc i;j � x k i;j respectively, and vehicle k stopping at node n or not can be determined by sn k n � sn o d i;n � v k i ð8i 2 P; n 2 NÞ. However, some rules must be complied with to ensure the feasibility of the route.

Rule 1:
Pd-pair j can connect to vehicle/pd-pair i when connect_to_judge i,j = 1.

Rule 2:
Each pd-pair transported by a vehicle must connect to the vehicle, or another pd-pair that connects after the third pd-pair or vehicle only once.

Rule 3:
Each vehicle/pd-pair must not be connected after by more than one pd-pair.

Rule 4:
Each pd-pair must not connect to itself.

Rule 5:
There cannot be circles in any route.

Rule 6:
Each pd-pair not being transported should not connect to any pd-pair or vehicle. The vehicle route must be as follows: pd-pair i 1 connects after vehicle k (namely x k pþ1;i 1 ¼ 1), pd-pair i 2 connects to pd-pair i 1 (namely x k i 1 ;i 2 ¼ 1), and pd-pair i 3 connects after pd-pair i 1 (namely x k i 1 ;i 3 ¼ 1), according Rule 2 and Rule 3. If x k i 2 ;i 3 ¼ 1, the length (e 5 +e 6 +e 7 ) between points d i2 and d i1 may be double-counted, because pd-pair i 2 does not connect after any pd-pair or vehicle. If x k pþ1;i 3 ¼ 1, the length (e 1 +e 2 +e 3 +e 4+ e 5 +e 6 +e 7 ) between point k and d i1 may be double-counted, as vehicle k has been connected after by pd-pair i 1 already.
All the rules have been considered in Section 3.3 (constraints 2, 3, 4, 5, 6). Fig 4 is a real-life connected graph, edges lengths are shown in it. In a feasible schedule, six pd-pairs (demand: 1) are transported by two vehicles (capacity: 5, distance limit: 10, stops limit: 6) along two routes. According the above definitions: Route 1 and route 2 are RSF; Pd-pair i 5 can be inserted into route 1 while pd-pair i 4 cannot; Pd-pair i 3 can be combined with i 1 while pd-pair i 4 cannot; Pd-pair i 2 and i 3 can connect to pd-pair i i , pd-pair i 4 can connect to pd-pair i 5 , and pd-pair i 5 can connect to pd-pair i 4 ; pd-pair i 3 can connect after pd-pair i 1 and i 2 , pd-pair i 4 can connect after pd-pair i 5 , each pd- Table 2. le e and lc i,j between two pd-pairs.

Instances
Order Paths Length of weighting sections (le e ) Length of connecting sections (lc i,j )

Fig 2(A)
j can connect to i p i -d i -p j -d j e 1 +e 2 +e 3 +e 7 +e 8 +e 9 L1 i cannot connect to j -e 7 +e 8 +e 9 +e 1 +e 2 +e 3 1

Fig 2(B)
j can connect to i p i -p j -d i -d j e 1 +e 2 +e 3 +e 4 +e 5 +e 6 +e 7 +e 8 +e 9 0 i cannot connect to j -e 4 +e 5 +e 6 +e 7 +e 8 +e 9 +e 1 +e 2 +e 3 1 The values of connect_to_judge i,j , connect_after_judge i,j and lc i,j for the routes combined with two pd-pairs are listed in S1 Appendix.  Table 3 (list non-zero variable only).

Neighborhood
Cordeau et al. [39] presented eight kinds of local search moves for OPDP: Couple-exchange, Block-exchange, Relocate-couple, Relocate-block, Multi-relocate, 2-opt-L, Double-bridge and Shake. Ribeiro et al. [40] and Ropke et al. [41] modified three large neighborhood removal heuristics and two large neighborhood insertion heuristics from Shaw P. [42], [43] and Potvin et al. [44] for OPDP. All those methods above may be of little efficiency for the OPDPSTRP in this paper, as the route structure of the OPDPSTRP is quite different from the classical OPDP, there are not so many pd-pairs can reinsert into a new route in the OPDPSTRP than in the classical OPDP. So five new neighborhood transform methods are presented for the OPDPSTRP.
Additionally, studies of Grimault et al. [45] and Ho et al. [7] show that solution feasibility of an OPDP is an important issue to the efficiency of the algorithm. So route structure feasibility judgement parameter of combining (RSFJPC, Definition 3 in Section 3.4) were applied in neighborhood transform methods (Insert, Spread, Point-delete) to ensure that each selected pd-pair can be inserted into the selected route. The evolution of RSFJPC will be studied in Section 4.3.

Insert.
In Insert, pd-pair i (may be carried by other routes, may not be carried, with pickup point p i and delivery point d i ) is selected randomly and inserted into a new route j chosen according pd_R_rs_judge i,j = 1.
As in Fig 5, pd-pair i is inserted into route j 2 from route j 1 , the scheme showed in Fig 5(A) is replaced by Fig 5(B). Stop point r 3 is removed from route j 1 because there is no pd-pair requiring transportation from/to r 3 , while structure of route j 2 remains the same. By reducing the number of stops, the route has been improved.

Spread.
In Spread, a pd-pair is selected and inserted into a new route as an Insert operation. Should the vehicle be overloaded, the success rate can be improved by choosing a new pd-pair i from route and transferring this into a new route j selected by pd_R_rs_judge i,j = 1, and this cycle will continue until the vehicle is no longer overloaded, or if the cyclic number k exceed the preset iterative numbers controlling value K. The task of preset value K is to control the computing time of this operation.
As illustrated in Fig 6(A), pd-pair i 1 transported through route j 1 is inserted into route j 2 , and the new scheme is shown in Fig 6(B). Stop point r 4 is deleted from route j 1 because there is no pd-pair requiring transportation from/to r 4 , and the structure of route j 2 remains the same. Since the vehicle is overloaded at section r 3 -r 4 in route j 2 , pd-pair i 2 is selected from route j 2 and inserted into j 3 , and then the new scheme is obtained as Fig 6(C), in which the vehicle is no longer overloaded at section r 3 -r 4 , while the structure of route j 3 remains unchanged, and the scheme is finally improved.

Point-delete.
Point-delete starts by choosing a route at random, before isolating the point with the minimal number picking stops and delivery stops on the route, and these pdpairs i2P are subsequently inserted into different routes j selected by pd_R_rs_judge i,j = 1, thus making it possible to delete the point from the first route.
As in Fig 7(A), Point-delete is operated on stop point r 4 in route j 1 , pd-pair i 1 and i 2 are inserted into route j 2 and j 3 respectively, and then the new scheme is obtained as Fig 7(B). Stop point r 4 is deleted from route j 1 because there is no pd-pair requiring transportation from/to it, while the structure of routes j 2 and route j 3 remains the same. The scheme is improved due to the reduction of stop point r 4 .

Route-delete.
In the Route-delete, net incomes ne i of each route i are computed; route k is selected according to the probability P k = ne k /∑ne i before being deleted from the scheme. All pd-pairs transported by route k are transferred to the state of non-carried. If there is no route with non positive net incomes, then the Point-delete should be executed.

Reassign-vehicle.
Reassign-vehicle is an Assignment Problem (AP) in which: rv i,j = 1 means route i being transported by vehicle j, and rv_benefit i,j is defined as its income. In this strategy, vehicles are reassigned to routes to achieve the best scheme by the Gurobi solver in Matlab.

Opt(k) and perturbation.
Since Reassign-vehicle may delivers the most significant change to the solution but cannot always improve the solution and requires more computer memory and time, it's better to be chosen according to a low probability to perturb the local best solution, so a kind of Perturbation is proposed to shock the local best solution instead of Reassign-vehicle. In Perturbation, Insert, Spread, Point-delete, Rout-delete and Reassign-vehicle are chosen according to the operators choosing probabilities p 1 , p 2 , p 3 , p 4 and p 5 respectively.

Route construction methods of neighborhood
Operators for routes construction and adjusting can be classed into three separate categories: pd-pair insertion, pd-pair deletion, and route deletion.

Route construction method for pd-pair insert.
Muelas et al. [28] arrayed two pdpairs by comparing the distance between pickup/delivery points of them based on four cases. Different from their methods, in this paper, a new approach is proposed to find the right locations in a new route and insert a pd-pair into it.
Firstly, it is evident that the section between any two neighbor stop points r i and r i+1 in route R is the shortest path since all pd-pairs must be transported along the shortest path. Once pd-pair k needs to be inserted into route R(r 1 −r 2 −r 3 . . .r n ), we have to find the right inserting location of pickup point p k and delivery point d k in route R first. Let l r i ;r j is the distance between r i and r j and let l Rðr i À r j À r k Þ ¼ l r i ;r j þ l r j ;r k is the length of section r i −r j −r k in route R.
A point r i meeting the condition l r i ;r iþ1 ¼ l r i ;p k þ l p k ;r iþ1 and a point r j meeting the condition l r j ;r jþ1 ¼ l r j ;d k þ l d k ;r jþ1 need to be found in route R. Finally, p k and d k need to be inserted into the Multi-start heuristic approaches for OPDP with shortest-path transport along real-life paths rear of the two points respectively, and we may get three types of results as follows: both r i and r j can be found, only one of r i and r j can be found, neither r i nor r j can be found.
Both r i and r j can be found: If both r i and r j can be found, and i<j, the structure form of the new route R' must be r−p−r −d−r. In order to ensure feasibility of the route structure, pd-pair k should be transported  along the shortest path in route R', that is, If both r i and r j can be found, and i = j, the route structure form of R' must be r−p−d−r. In order to ensure feasibility of the route structure, r i −p k −d k must be the shortest path, that is, Otherwise, pd-pair k cannot be inserted into route R. As in Fig 9, l R 0 ðr i À p k À d k Þ 6 ¼ l r i ;d k , so the route structure is unfeasible.

➢ Other cases
If both r i and r j are found, and i>j, pd-pair k are opposite to route R and cannot be inserted into it, as in Fig 10. Because there is only one shortest path between any two nodes in this paper, if there is more than one r i or r j satisfying the criteria, it must be sure that p k or d k is in route R already, then point p k or d k need not be inserted into it repeatedly. As in Fig 11, l r Only one of r i and r j can be found: If only r i can be found, the route structure form of R' must be r−p−r−d. Pd-pair k must be transported along the shortest path in route R', that is, l R 0 ðp k À r iþ1 ...r iþ3 À d k Þ ¼ l p k ;d k . Otherwise, pdpair k cannot be inserted into R. As in Fig 12, l R 0 ðp k À r iþ1 ...r iþ3 À d k Þ ¼ l p k ;d k , so the route structure is unfeasible.

➢p−r−d−r
If only r j can be found, the route structure form of R' must be p−r−d−r. Pd-pair k must be transported along the shortest path in route R', that is, l R 0 ðp k À r jÀ 2 ...r j À d k Þ ¼ l p k ;d k . Otherwise pdpair k cannot be inserted into route R. As in Fig 13, l R 0 ðp k À r jÀ 2 ...r j À d k Þ 6 ¼ l p k ;d k , so the route structure is unfeasible.
Neither r i nor r j can be found: If none of r i and r j can be found, and l R 0 ðpÀ r 1 ...r n À dÞ ¼ l p k ;d k , the route structure form of route R' must be p−r−r−d, pd-pair k is transported along the shortest path in route R', and the route structure must be feasible, that is, pd-pair k can be inserted into route R. As in Fig 14, l R 0 ðpÀ r 1 ...r n À dÞ 6 ¼ l p k ;d k , so the route structure is unfeasible.
If neither r i nor r j are found, and l R 0 ðpÀ r 1 ...r n À dÞ 6 ¼ l p k ;d k \ l d k ;r 1 <¼ l r n ;p k , pd-pair k must be connected to the front of route R, the route structure form of R' must be p−d−r−r to minimize the length of R'. Set path(R) is the vehicle path which consists of all pass points in route R. Pd-pair k must be transported along the shortest path and each point must not be visited more than once in route R' to ensure feasibility of the route structure, that is, Otherwise pd-pair k cannot be inserted into route R. As in Fig 15(A), the route structure is unfeasible because l p k ;r 1 6 ¼ l p k ;d k þ l d k ;r 1 . As in Fig 15(B), the route structure is unfeasible because If neither r i nor r j are found, and l R 0 ðpÀ r 1 ...r n À dÞ 6 ¼ l p k ;d k \ l d k ;r 1 > l r n ;p k , pd-pair k must be connected to the back of route R, and the route structure form of route R' must be r−r−p−d to minimize the length of route R'. Pd-pair k must be transported along the shortest path and each point must not be visited more than once in route R' to ensure feasibility of the route structure, that is, l r n ;d k ¼ l r n ;p k þ l p k ;d k and l pathðR 0 nÀ 1 Þ;p k ¼ l pathðR 0 nÀ 1 Þ;pathðR 0 n Þ þ l pathðR 0 n Þ;d k . As in Fig 16  (A), the route structure is unfeasible because l r n ;d k 6 ¼ l r n ;p k þ l p k ;d k . As in Fig 16(B), the route structure is unfeasible because l pathðR 0 nÀ 1 Þ;p k 6 ¼ l pathðR 0 nÀ 1 Þ;pathðR 0 n Þ þ l pathðR 0 n Þ;d k . 4.2.2 Route construction method for pd-pair delete. Point p k or d k needs to be deleted if there is no other pd-pair that needs picking up from it or delivering to it after deleting pd-pair k from route R. Otherwise, it should remain in route R.

Route construction method for route delete.
For Route-delete strategy, all pd-pairs are removed and route is deleted. quickly, some theories will be studied as follows.
Theorem 1: A necessary and sufficient condition of which pd-pair k can be inserted into route m according route structure is that pd-pair k can be combined with all pd-pair transported by route m according route structure.
Lemma 1: If pd-pair i is inserted into route m, and a new route m' is obtained, pd-pair j which cannot be combined with pd-pair i according route structure cannot be inserted into route m' according route structure. Multi-start heuristic approaches for OPDP with shortest-path transport along real-life paths Lemma 2: Assuming pd-pair k cannot be inserted into route m according route structure. If all pd-pairs which cannot be combined with pd-pair k according route structure are deleted from route m, then pd-pair k can be inserted into the new route m' evolving from route m according route structure.
Lemma 3: A new route R m can be constructed with pd-pairs according route structure, which can be combined with each other according route structure. The proofs of Theorem 1, Lemma 1, Lemma 2, and Lemma 3 are provided in S2 Appendix.

Evolution methods for route structure feasibility judgement matrix.
As RSFJMC will remain the same at each local search move step, while RSFJMI will change as routes evolve. If RSFJMI is updated at each step, too much computing time will be spent on carried out the algorithm. In this section, evolution methods of RSFJMI will be studied to ensure that the calculation time is acceptable.
According to the route construction methods in Section 4.2, the evolution of route structures can be grouped into two categories: pd-pair deletion and pd-pair insertion. The RSFJMC and RSFJMJ can be calculated after an initial feasible solution is generated. At each local search move step, the RSFJMI will be updated by methods as follows.

Evolution of RSFJMI for pd-pair insertion.
When pd-pair i is inserted into route k, and a new route k' is acquired. In order to reduce calculating time, let L 1 = {l|[pd_combi-ne_rs_judge](l,i)==0}, [pd_R_rs_judge](:,k') for route k' can be updated by formula (18) according to Lemma 1.

Evolution of RSFJMI for pd-pair deletion.
When pd-pair i is deleted form route k, and a new k' is acquired. In order to reduce calculating time, Let L 1 = {l|[pd_combine_rs_judge](l,i)==0}, and let L 2 = {p l −d l }, a set of pd-pairs transported by route k'.
Additionally, Lemma 3 can be used when new routes are generated, such as for the generation of new routes.

Generation of initial solution
It is important to obtain a higher-level initial solution for such a large-scale and complex problem. In this paper, a generation method of initial solution is proposed, based on the idea of maximum saving. p i and d i are defined as pickup point and delivery point of pd-pair i correspondingly. l(k,l) is defined as the length between point k and point l.
Generation steps of initial solution are presented in Algorithm 1 (Fig 17).

Algorithm
Many Local Search (LS) meta-heuristics are studied for the VRP, such as the ALNS proposed by Ropke [59]. The VND, a kind of local search in which wider and wider neighborhoods are successively explored, is studied by Salehipour et al. [60]. Variable Neighborhood Search (VNS), which is similar to the VND but with perturbation, is studied by Hansen et al. [61], which can reduces the impact of the initial solution on the algorithm performance (Brimberg et al. [62]). Since the major neighbourhood operator is pd-pair insertion and a choosing pd-pair need be inserted a fixed position in a route at each local search move step in the OPDPSTRP, solutions are not easy to be changed and always cannot be improved by only one operator obviously. That is to say that new methods which is different from traditional LS need be proposed for the OPDPSTRP. A basic VND and a basic VNS are proposed for the OPDPSTRP based on the above five neighborhoods. A new MS_VND and a new MS_VNS are developed to improve the efficiency of the proposing neighborhoods and algorithms in this paper.

VND and VNS.
A VND and a VNS are processed to solve the OPDPSTRP, Pseudocode of them are presented as in Algorithm 2 (Fig 18) and Algorithm 3 (Fig 19).

MS_VND and MS_VNS.
In order to achieve a near optimal solution of this problem, we propose a new MS_VND metaheuristic and a new Multi-Start Variable Neighborhood Search (MS_VNS) metaheuristic.
The steps of MS_VND in this paper are presented as in Algorithm 4 (Fig 20).
The steps of MS_VNS in this paper are presented as in Algorithm 5 (Fig 21). As in Algorithm 4 and Algorithm 5, the MS_VND and the MS_VNS have been improved in the following ways: (1) A Multi-start candidate solution set with size of n is acquired from the initial solution and updated to diversify the search, and parts of worse candidate solutions are replaced by the best solution according replacing proportion m if the solution is improving. (2) Five new operators (opt1, opt2, opt3, opt4, opt5) are utilized to improve the solution, which is different from traditional VNS. (3) In MS_VND, each candidate solution is transformed only once at a step (Multi-start-candidate-solution and One-operator), which is different from traditional LS (One-candidate-solution and Local search). (4) A new local solution inferior to the primary one can also be accepted on the basis of two hypotheses: the iterations Multi-start heuristic approaches for OPDP with shortest-path transport along real-life paths of that the solution keep the same are more than half of the presetting value, and the evaluation of the new solution is not drastically worse than the primary one.
Changes and complexity ranking (from low to high) of operators can be primarily concluded as: Insert<Spread<Point-delete<Route-delete<Perturbation according to the features of them. And this ranking will be tested further in Section 5.

Instances and computational results
To evaluate the performance of the VNS, the VND, the MS_VND and the MS_VNS, we designed 84 of the OPDPSTRP instances partitioned in small size connected graphs (3×4), medium size connected graphs (6×8) and large size connected graphs (10×10) at random. The Gurobi MIP solver (version 7.5.2) was used to compare with these approaches proposed in this paper. The following sections will describe the generation of instances, list the parameter values used in these approaches, and provide test results and optimality gaps.

Generation of instances
Each instance name has a format m0×m1−m2−m3−m4−L/H, in which m0×m1 is the size of a connected graph, each edge in this graph is deleted according to 1/m2, 1/m3 is the pd-pairs generation probability between each two nodes, 1/m4 is the vehicle generation percentage for each node and L/H is low or high cost. Consider instance 3×4-10-3-3-L as an example: the size of the incomplete digraph is 3×4(12 nodes and 144 node-pairs), the deletion probability of each edge is 1/10, pd-pair generation probability between each two nodes is 1/3, the vehicle generation percentage for each node is 1/3 and L means low cost. It has been checked that there is only one shortest path between any two nodes in each graph. Data of these 84 instances can be found in S3 Appendix.

Parameter setting
All experiments were conducted on a desktop equipped with an Intel(R) Core(TM) i7-4510U 2.00 gigahertz processor and 8 Gigabyte RAM. The operating system of this PC is a 64-bit Window 8. The Gurobi solver 7.5.2 was used to solve the MIP model, and all algorithms in this paper was coded in Matlab R2015a.
The MIP model is solved by the Gurobi solver with termination conditions wherein computing time is over 108000 seconds or the Gurobi_Gap (the gap between the best feasible objective value and the upper bound of the MIP model, which will be introduced in Section 5.3) is less than 5%. The long preset time aims to ensure that the Gurobi solver can obtain at least one feasible solution served as a comparison indicator with the proposed approaches, although in some cases it failed to achieve this goal. We have provided the upper bounds found by Gurobi solver as well, for a more in depth reference to the performance of the proposed approaches.
The parameter setting of the algorithms have been tested based on 9 instances mention in S4 Appendix (including small size connected graphs, medium size connected graphs and large size connected graphs), the results are studied in S4 Appendix. The algorithms parameter setting is tuned by determining a trade off between solution quality and CPU time after numerous experiments.
According to S4 Appendix, the operators sequence is determined as Insert, Spread, Pointdelete, Rout-delete and Perturbation (k = 1, 2, 3, 4 and 5) in these approaches, and the values of the parameters are gathered in Table 4.

Test results
The abbreviations of the experimental indicators and corresponding definitions are listed in Table 5.
Results solved by the Gurobi solver, the VND, the VNS, the MS_VND and the MS_VNS are shown in Table 6 (small size graphs), Table 7 (medium size graphs) and Table 8 (large size graphs). Each instance has been solved 10 times by each algorithm.
Comparison of the average calculation efficiency of the VND, the VNS, the MS_VND and the MS_VNS are shown in Figs 22-24.
The above results shows that: 1. CPU time of the Gurobi solver, the VND, the VNS, the MS_VND and the MS_VNS are independent of road network scale. Numbers of pd-pairs and vehicles are the major influential factors to the CPU time of the Gurobi solver, numbers of pd-pairs are the major influential factors to the CPU time of the VND, the VNS, the MS_VND and the MS_VNS.
2. In general, instances with no more than 100 pd-pairs can get an optimum solution by the Gurobi solver in preset termination CPU time (108000 seconds). Instances with more than 5000 (pd-pairs×vehicles) always cannot receive the optimal solution with no more than 20% Gap_Gurobi within an acceptable CPU time (108000 seconds) by Gurobi solver. Large instances with over 40000 (pd-pairs×vehicles) always cannot get feasible solutions within 108,000 seconds by the Gurobi solver. 3. In almost all instances, the VND, the VNS, the MS_VND and the MS_VNS with the operators proposed in this paper always can acquire the optimal solution with no more than 10% Gap to the Gurobi solver within an acceptable CPU time. The MS_VND can acquire the optimal solution with no more than 5% Gap to the Gurobi solver within an acceptable CPU time. The MS_VNS always takes a little more CPU time than the MS_VND to solve the instances with more pd-pairs. 4. The VND, the VNS, the MS_VND and the MS_VNS can even get better solutions than the Gurobi solver for the instances with large numbers of pd-pairs and vehicles.
5. In almost all instances, the MS_VND significantly outperforms the VND, the VNS and the MS_VNS in terms of average solution quality (Figs 22-24). The MS_VND acquires almost all the best average solutions and most of the best solution for each instance (boldface numerical value in Table 6, Table 7 and Table 8).

Conclusions
A new Pickup and Delivery Problem with new route structure, the OPDPSTRP, which is proposed in real-life connected graphs, is introduced and formulated in a new way. Five operators Table 5. Abbreviation of experiment indicators and definitions.

Gurobi_UB
The upper bound of the MIP model obtained by the Gurobi solver in a preset running time.

Gurobi_LB
The best feasible objective value found by the Gurobi solver in a preset running time.

VND_LB_Avg
The average feasible objective value obtained by the VND after a preset number of iterations.

VNS_LB_Avg
The average feasible objective value obtained by the VNS after a preset number of iterations.

MS_VND_LB_Avg
The average feasible objective value obtained by the MS_VND after a preset number of iterations.

MS_VNS_LB_Avg
The average feasible objective value obtained by the MS_VNS after a preset number of iterations.

VND_LB_Best
The best feasible objective value obtained by the VND after a preset number of iterations.

VNS_LB_Best
The best feasible objective value obtained by the VNS after a preset number of iterations.

MS_VND_LB_Best
The best feasible objective value obtained by the MS_VND after a preset number of iterations.

MS_VNS_LB_Best
The best feasible objective value obtained by the MS_VNS after a preset number of iterations.

VND_Time
Average CPU time for solving the MIP model by the VND (second).

VNS_Time
Average CPU time for solving the MIP model by the VNS (second).

MS_VND_Time
Average CPU time for solving the MIP model by the MS_VND (second).

MS_VNS_Time
Average CPU time for solving the MIP model by the MS_VNS (second).     Table 5 https://doi.org/10.1371/journal.pone.0227702.t006 Table 7. Computational results for instances of medium size graphs.      Table 5. https://doi.org/10.1371/journal.pone.0227702.t008 Multi-start heuristic approaches for OPDP with shortest-path transport along real-life paths are proposed for the OPDPSTRP. A basic VND, a basic VND, a new MS_VND, and a new MS_VNS are developed for this problems and compared with the Gurobi solver. The test results show that the VND, the VNS, the MS_VND and the MS_VNS always can acquire the optimal solution with no more than 10% Gap to the Gurobi solver within an acceptable CPU time. In almost all instances, the MS_VND significantly outperforms the VND the VNS and the MS_VNS in terms of solution quality. The MS_VNS always takes a little more CPU time than the MS_VND to solve the instances with more pd-pairs. In the large instances, the MS_VND is still able to generate good feasible solutions in a reasonable CPU time, which is of vital practical significance for real-life instances.

Instances
Our future work will concentrate on three aspects: (1) In order to improve the algorithm efficiency of the OPDPSTRP, path feasible strategy for neighborhoods will be studied. (2) There is only one shortest path between any two nodes in the connected graph in this paper, but there may be more than one shortest path between two nodes in some other real-life road networks, which we are going to study in the future. (3) The OPDPSTRP with time windows will be studied in the future.