To perform the optimisation-based analysis, we need to set up suitable optimisation models and solution methods. An integer linear programming model (ILP) used here takes the specific cost structure of our real world case into account. It is an adaption of a model presented by Desaulnier et al. (2010) adding a distinction between the two different vehicle types. To tackle larger instance, also a meta-heuristic solution approach is proposed.
Integer linear programming model (ILP). We tackle the optimisation problem by formulating a mathematical model that can be solved using integer linear programming solvers. At the same time, this provides a formal problem description. For are given a set of nodes V=B∪’P ∪’D where B denotes depot nodes, P denotes pickup nodes, and D denotes delivery nodes the complete graph G=(V,E) describes connections in an underlying road network. The travel time between two nodes i,j∈V is given by ti,j∈R. For a pickup (delivery) node i∈V its corresponding delivery (pickup) node is denoted by vi∈V. The fleet of vehicles K=Ks∪’Kf is composed of standby vehicles Ks and flexible vehicles Kf, which differ in their usage costs. For each vehicle k∈K, we are given its start location sk∈B and end location ek∈B, and its capacity Ck∈N. Costs differ for loaded and unloaded vehicles, so costs per distance for loaded vehicles are denoted by csl∈R+ and cfl∈R+, and for empty vehicles by cse∈R+ and cfe∈R+, with cse≤csl and cfe≤cfl. For each node i∈V, we are given a time window [ai,bi] and a service duration si∈R+ reflecting the desired departure time of passengers and a travel duration which include a buffer time that is varied within the numerical experiments. For each node i∈V, its load li∈N specifies the number of passengers that want to get in (li>0) or out (li<0) at i. A sufficiently large constant M∈K is introduced for modelling purposes. The following integer linear programming model is an adaption of the formulation given in (Desaulnier et al. 2010) to the problem at hand. It uses the following variables. Binary variables xik,j describe if an arc (i,j∈E) is used by the vehicle k∈K. Binary variables xik,jl specify if an arc (i,j∈E) is used by a loaded vehicle k∈K. Binary variables xik,je specify if an arc (i,j∈E) is used by an empty vehicle k∈K. Variables Ti∈R determine for nodes i∈N the start time of the service. Variables Li∈R determine for nodes i∈N the load (in the number of passengers) when a vehicle leaves that node.
The formulation of the model is provided in the figure below. The objective function (1) reflects the aim of minimising costs, which depend on the load and type of the vehicle. Constraints (2) to (11) are taken from (Desaulnier et al. 2010) for modelling classical dial-a-ride constraints. In addition, constraints (12) and (13) are introduced for distinguishing loaded from empty vehicles. The general purpose mixed integer linear programming solver IBM ILOG CPLEX Optimizer, version 12.6.2, was used. For resolving larger instances, a metaheuristic approach based on a best insertion construction heuristic and a subsequent variable neighbourhood descent was applied.

Variable Neighborhood Search (VNS). With the integer linear programming formulation, we are only able to solve instances up to 36 requests within the time limit of two hours. For tackling larger instances, we developed a metaheuristic approach using the framework of a general variable neighbourhood search (VNS) (Mladenović & Hansen, 1997). The goal of the VNS is to compute (near-) optimal routes in a short amount of time, however, without giving any guarantees about optimality. In the developed VNS, we first constructed an initial solution by following the best insertion strategy, which starts from a set of empty routes and iteratively adds a request (in an arbitrary order) to the route and position that result in the least cost increase. The resulting solution is then further improved via a variable neighbourhood descent (VND) algorithm using a nested neighbourhood structure of cardinality four. These four considered neighbourhood structures, each defining types of allowed moves, i.e., small changes from the starting solution, are intra-relocate, exchange, inter-relocate, and block-relocate. The relocate neighbourhood structure re-assigns a request by moving it from one position on the first trip to another position on a second trip. These positions can be on the same trip (intra-relocate), or within different tours (inter-relocate), yielding two variants of that neighbourhood. In the exchange neighbourhood structure, the positions of two requests from two different trips are swapped. A move in the block-relocate neighbourhood changes the position of a consecutive sequence of requests. After no improvement can be found in these neighbourhoods anymore and a local optimum is reached, a shaking procedure makes random changes to the incumbent solution and the VND starts anew. The magnitude of these changes increases for each iteration of the VND, in which no improvement was found after eight moves in the inter-relocate neighbourhood. The VNS terminates after a given number of iterations or a given amount of time, yielding the best solution found during the search.
Assessment of Solution Methods. To determine the quality of the solution methods, we conducted a range of numerical experiments. Table 1 provides an overview. We see that the integer linear programming approach is only suitable for small instances. The gap of the VNS approach to the best solutions found is very small and well below one per cent on average. Though the VNS is a heuristic approach, the results of Table 1 confirm that implemented VNS approach yields high-quality solutions. The figure below visualises two examples of solutions obtained by the two implemented solution methods.
Table 1. Comparison of Variable Neighborhood Search with results of the ILP model.
#locations | percentage | #solved exactly | VNS gap to best |
<20 | 7.0% | 80.6% | 0.30% |
20-39 | 15.6% | 65.6% | 0.32% |
40-59 | 3.5% | 38.9% | 0.23% |
60-79 | 8.2% | 11.9% | 0.08% |
80-99 | 11.7% | 0.0% | 0.00% |
100-119 | 17.1% | 0.0% | 0.00% |
120-139 | 24.5% | 0.0% | 0.00% |
140-159 | 10.1% | 0.0% | 0.00% |
160-180 | 2.3% | 0.0% | 0.00% |

Next: Cross-border transfer of a systemic Mobility-as-a-Service -related innovation
References
G Desaulniers, J Desrosiers, A Erdmann, MM Solomon, and F Soumis. VRP with Pickup and Delivery. In The vehicle routing problem, pages 225-242. Society for Industrial and Applied Mathematics, 2001.
Mladenović, N., and P. Hansen. Variable neighborhood search. In Computers & Operations Research 24 (11), pages 1097-1100, 1997.