Solving a Large-Scale Multi-Depot Vehicle Scheduling Problem in Urban Bus Systems

This study proposes an improved model and algorithm for the large-scale multi-depot vehicle scheduling problem (MDVSP) with departure-duration restrictions. In this study, the time-space network is applied to model the large-scale MDVSP. Considering that crews usually change shifts in the depot, departure-duration restrictions are added to the classic set-partitioning model to ensure that buses return to the depot when crews reach their working time limits. By embedding a preliminary exploring tactic to the shortest path faster algorithm (SPFA), researchers developed an improved large neighborhood search (LNS) algorithm to solve large-scale instances of MDVSPwith departure-duration restrictions.The proposedmethodology is applied to a real-life case in China and several test instances. The results show that the improved LNS algorithm can achieve very good performance in computational efficiency without deteriorating solution quality, which is important for large-scale systems. More specifically, the total cost of the improved LNS algorithm is approximately equal to branch-and-price, but the computational time is much shorter in the case study. For test instances with different number of timetabled trips (500, 1000, 1500, and 2000), the Quality Gap (QG) is very small, approximately 0.35%, 0.38%, 0.63%, and 0.93%, while the Efficiency Ratio (ER) reaches up to 2.89, 2.98, 3.65, and 3.79, respectively.


Introduction
Bus operation planning commonly includes four components: (1) bus network design; (2) timetable design; (3) bus scheduling; and (4) crew scheduling [1].Among them, bus scheduling is extremely complex.An efficient scheduling plan is of great benefit to both bus companies and passengers.Compared to the single-line scheduling mode universally used in China, the multi-depot vehicle scheduling mode which belongs to the regional scheduling mode is more efficient.In the multi-depot vehicle scheduling problem (MDVSP), buses are allocated globally, helping bus companies reduce operating costs and promote operational efficiency.Although several models and algorithms have been developed for the vehicle routing problem (VRP) or vehicle scheduling problem (VSP), the extremely large-scale MDVSP involving thousands of timetabled trips and several depots remains relatively unexplored.Most algorithms are not efficient enough to provide a relatively optimal solution within an acceptable computational time.Moreover, bus scheduling systems in China are particularly large and busy due to the large population base and city size.In 2013, there were 21,293 vehicles running on 754 bus lines in Beijing every day, and, by 2014, the number of bus lines had increased to 652 in a third-tier city like Ningbo.For these large-scale systems, it is considerably difficult to apply the multi-depot vehicle scheduling mode into practice.The key issue of application lies in the availability of methods to downsize the scheduling system and improve the algorithm.As background, previous studies are analyzed from three aspects (i.e., network description, model formulation, and solution algorithm) as follows.
According to previous studies, both connection-based networks and time-space networks could be applied to model the VSP.In early studies, vehicle scheduling problems were modeled as connection-based networks [2].The time-space network was first proposed for air scheduling due to its simplicity in modeling possible connections between flights [3].Kliewer et al. were the first to apply the time-space network to VSP; they then proposed a commodity network flow model on the basis of the time-space network structure [4,5].For a small-scale VSP, the connection-based network was superior, whereas, for a middle or large-scale VSP, the time-space network could achieve better performance in downsizing the system by reducing the number of deadhead arcs [6].Gintner et al. and Naumann et al. efficiently solved a practical large-scale and middle-scale MDVSP, respectively, using the commodity network flow model based on the timespace network [7,8].However, the reduction of deadhead arcs in the time-space network has not been sufficiently considered in previous studies.Instead of checking and reducing deadhead arcs after the network takes its shape, we avoid generating unnecessary deadhead arcs in the network building process to achieve the goal of reducing the number of deadhead arcs.In this way, we can reduce not only the number of deadhead arcs but also the time used in constructing the time-space network so as to improve efficiency.
For the VRP, there are three basic models from previous studies (i.e., vehicle flow model, commodity network flow model, and set-partitioning model) and several variations involving different constraints (e.g., VRP with time windows [9], VRP with backhauls [10], and VRP with pickup and delivery [11]).The goals of these models normally fall into two categories, minimizing fleet size or minimizing operating costs.The set-partitioning model was originally proposed by Balinski and Quandt [12].As the route feasibility is implicitly considered in the definition of the circuit set, the set-partitioning model can easily take extra constraints into account [13].In recent years, there have been too many safety incidents related to bus drivers working overtime in China.Drivers working overtime not only jeopardize their own health but also threaten the safety of passengers.Therefore, compared to other constraints, crews' working time limits are the primary consideration taken into account in this study.Drivers should strictly adhere to the change-of-shift time in order to get sufficient rest, which will ensure the safety of passengers on board.As drivers usually change shifts in the depot, buses should return to the depot when the crews reach their working time limits to ensure safe driving.Based on these conditions, the MDVSP in this study is formulated as a set-partitioning model with departure-duration restrictions.
As VRP is an extension of the Traveling Salesman Problem (TSP) and the VSP explored in this study is a VRP that arises in the public transport area, the algorithms for VSP can build on the extensive and successful work done for TSP and VRP.The algorithms in previous studies for VRP are categorized into two groups: exact algorithms and heuristic algorithms.Exact algorithms can achieve high quality solutions when dealing with small-scale systems; however, they suffer great limitations in practice.In contrast, the heuristic algorithms that usually perform a relatively limited exploration of the search space are capable of solving larger VRP problems within modest computational time.However, for extremely large-scale real-life scheduling instances involving several thousand trips, most of the available algorithms are computationally intensive, requiring thousands of CPU seconds.The branch-and-bound algorithm belongs to the exact algorithms category and has been extensively applied to solve the VRP and its variants in recent decades [14][15][16].As noted by Toth and Michalewicz, branch-and-bound could find a significantly high quality solution but was not capable of solving large instances encountered in practice [13,17].The column generation was an important technique that could solve larger Linear Programing (LP) problems [18].Many scholars have applied column generation to solve VRP [19,20].By integrating column generation techniques with the branch-and-bound scheme (which is also known as branch-and-price), many systems can be solved with significantly good quality solutions, as demonstrated by many successful studies in recent years [21][22][23][24][25].However, the computational efficiency of the branch-and-price is not so satisfying, in particular for large-scale systems.As the MDVSP explored in this study is extremely large, a heuristic algorithm is established on the large neighborhood search (LNS) framework.By integrating branch-and-price into LNS, an algorithm characterized by both high quality solutions and computational efficiency is established.As the MDVSP in this study has departure-duration restrictions, the shortest path faster algorithm [26] is modified by embedding a preliminary exploring tactic to solve the subproblem.
In contrast to previous studies on MDVSP, the contributions of this study are as follows: (1) presenting a novel way to deal with deadhead arcs so as to improve the efficiency of downsizing the time-space network based MDVSP, (2) modeling the MDVSP by adding departureduration restrictions to ensure safe driving which is imperative in China, (3) proposing an improved LNS algorithm to ensure good performance both in solution quality and computational efficiency, and (4) modifying the shortest path faster algorithm (SPFA) by embedding a preliminary exploring tactic to address departure-duration restrictions.
The remaining sections are organized as follows.Section 2 introduces basic descriptions of the MDVSP and network modeling.Sections 3 and 4 present the model formulation and solution algorithm to solve a large-scale MDVSP with departure-duration restrictions.Sections 5 and 6 evaluate the performance of the proposed methodology using a real-life case in China and several test instances.The conclusions and perspectives are summarized in the last section.

Problem Statement and Network Modeling
Before constructing the mathematical formulation, the timespace network is first introduced.In this section, basic descriptions and assumptions are presented.The methods of network modeling are explored at length.Some necessary definitions and terms are defined in this section.

Basic Descriptions and Assumptions.
The VSP is a VRP encountered by bus companies when addressing the task of assigning buses to cover a given set of timetabled trips with the goal of minimizing fleet size or operating costs.Each timetabled trip is accomplished exactly once by a bus, and each bus performs a feasible sequence of trips.The MDVSP in this study involves several depots, and a timetabled trip can be run by any bus from any depot.The large-scale MDVSP with departure-duration restrictions as addressed in this study is built upon the following hypotheses: (1) The scheduling scheme is organized in terms of individual days, and minutes are the smallest unit of time.
(2) All buses are homogeneous and the number of buses is limited.
(3) Every bus belongs to a home depot.At the end of the day, a bus has to return to the same home depot where it started in the morning.
(4) The cost for buses waiting outside the depot is very high.It is more favorable to wait at a depot than at other stations.If there is enough time, the bus should return to the depot and wait.

Network Description.
The VSP can be modeled by either a connection-based network or a time-space network.Using the connection-based network, the problem can only be slightly downsized as the reducible long arcs are in a relatively small quantity.In contrast, there are significant numbers of deadhead arcs that can be removed in the time-space network.If the problem contains m terminal stations and n trips, the number of deadhead arcs in a time-space network can be reduced to o(nm), in comparison to o(n 2 ) in the connectionbased network.In general, for a large-scale problem, the number of timetabled trips n is much larger than the number of terminal stations m.Consequently, the problem scale is much smaller when using a time-space network compared to using a connection-based network for the same scheduling system.Moreover, as each timetabled trip is characterized by a departure time and arrival time as well as an origin and destination station (time and space attributes), the time axis and space axis in the time-space network constitute a two-dimensional space which can more clearly identify the two attributes of each timetabled trip and the time-space relationships of different bus tours.On that account, it is more favorable to use a time-space network to handle the largescale MDVSP as done in this study.Figure 1 illustrates the time-space network of depot k in MDVSP.
The time-space network is a directed graph composed of many vertices and arcs.All vertices have two attributes (time and space) and each of them connects a group of possible arrival arcs to the possible departure arcs of the following group.Each arc corresponds to a transition in time or space.Every depot and terminal station vertex indicates a possible arrival or departure event (an arc) at the depot or a given station at a certain time.For instance, s11 denotes a bus that is reaching and then leaving terminal station s1 to run a timetabled trip (s11, s21) at a time around 6:00.
For a depot k ( ∈ , K is the set of all depots), let   =(  ,   ) be a complete graph, where   = {(), ()} ∪  ∪  is the vertex set.o(k) corresponds to the set of depot vertices where buses pull out from depot k at a certain time(e.g., k1, k2, and k3 in Figure 1), whereas d(k) corresponds to the set of depot vertices where buses pull into depot k at a certain time (e.g., k5, k8, and k9 in Figure 1).T is the set of terminal station vertices where tasks start at a certain station and time (e.g., s11, s12, and s13 in Figure 1), whereas E is the set of terminal station vertices where tasks end at a certain station and time (e.g., s21, s22, and s23 in Figure 1).There are five types of arcs in arc set   , including pull-in arcs, pull-out arcs, deadhead arcs, waiting arcs, and task arcs (as shown in Figure 1).A pull-out arc (o,  ),  ∈ (),   ∈ , denotes a bus pulling out from depot k in order to start a sequence of timetabled trips.A pull-in arc (  , d),   ∈ ,  ∈ (), denotes a bus pulling into depot k after completing a sequence of timetabled trips.A task arc (  ,   ),   ∈ ,   ∈ , represents a timetabled trip.After serving a timetabled trip, each bus can wait to serve one of the trips starting later from the same station where it is holding, or it can change its location by moving unloaded to another station in order to serve the next loaded trip starting there.The former is a waiting arc, which indicates a bus that waits at a place for some time.There are two kinds of waiting arcs depending on whether the bus waits at a terminal station or at the depot.The latter is an unloaded trip corresponding to the deadhead arc.For a deadhead arc (  ,   ),   ∈ ,   ∈ ,   is an end vertex in set E and   is a start vertex in set T. For a pair of task arcs (1  , 1  ) and (2  , 2  ) that end at station vertex 1  at time t a and start at station vertex 2  at time t b , respectively, if   +   <   , where t ab is the travel time of a bus from station 1  to 2  , task arc (1  , 1  ) and (2  , 2  ) are termed as a compatible pair of tasks.The compatibility relationship represents whether trip (2  , 2  ) can be served after trip (1  , 1  ) by the same bus.Any pair of compatible task arcs can be connected by a deadhead arc showing a possible tour of a bus.For example, a pair of compatible task arcs (s11,s21) and (s42,s32) are connected by a deadhead arc (s21,s42) in Figure 1, and the deadhead arc (s21,s42) indicates that a bus moves unloaded from station s2 to station s4 at time t s21 and then waits at station s4 until t s42 .Figure 1 is the time-space network of depot k.For multiple depots (K), there are K layers of Figure 1 that stack up together with shared timetabled trips (task arcs).
Most deadhead arcs can be removed from the timespace network.Among all the tasks compatible with task arc (1  , 1  ), the nearest one is designated as the first match arc.Let D be the set of deadhead arcs that connect task arc (1  , 1  ) with all its compatible task arcs.Deadhead arc f,  ∈ , is the connection between task arc (1  , 1  ) and its first match arc.Other arcs apart from f in set D can be removed as each of them can be represented by the deadhead arc f and a waiting arc.The total number of arcs will be reduced significantly compared to the original situation.Nevertheless, all possible connections remain feasible.For example, task arcs (s42,s32) and (s43,s33) are both compatible with task arc (s11,s21) in Figure 1, whereas task arc (s42,s32) is the first match arc with arc (s11,s21).Deadhead arc (s21, s43) can be replaced by deadhead arc (s21, s42) and waiting arc (s42, s43).Let T be the set of task arcs that end at terminal station s1 with the same first match task arc t that starts at terminal station s2.If task arc t1∈T is the nearest task arc to t among all elements in T, t1 is designated as the nearest match arc of t.For instance, in Figure 2, arc t4 is the first match of arc t1, t2, and t3.As t3 is nearest to t4; t3 is the nearest match arc of t4.Deadhead arc (s11,s21) in Figure 2 can be removed as it can be replaced by deadhead arc (s13,s21) and waiting arc (s11,s13).Similarly, deadhead arc (s12,s21) can be replaced by deadhead arc (s13,s21) and waiting arc (s12,s13).Only deadhead arc (s13,s21) is retained.
Each arc is associated with a nonnegative cost.
where   ,   ,    ,   ,   , and   in Figure 1 denote the cost of a pull-out arc, pull-in arc, waiting arc outside the depot, waiting arc inside the depot, task arc, and deadhead arc, respectively; ,   ,   , and  represent the unit cost of a bus running unloaded and loaded and waiting outside the depot and inside the depot, respectively; l ds , l sd , and l ss denote the distance from the depot to a station, from a station to the depot, and from a station to another station, respectively; t s , t d , and t t represent the waiting time at a station and at the depot and the duration time of a timetabled trip.A deadhead arc cost may include the waiting cost as a bus may arrive early before the next timetabled trip starts.
In the time-space network, the MDVSP can be described as many circuits (bus tours) starting from their respective depot vertex, passing by different kinds of arcs, and eventually ending at their depot vertex.Every task arc is covered exactly once.Every circuit, consisting of a feasible sequence of different types of arcs chained with each other, starts with a pull-out arc and ends up with a pull-in arc, which means a bus leaves the depot to start performing timetabled trips and goes back to the depot after finishing all timetabled trips.Take a circuit "k1-s11-s21-k5-k6-k7-s43-s33-k12" in Figure 1, for example: the bus first drives out from depot  (pull-out arc (k1, s11) to station s1 to perform a timetabled trip (task arc (s11, s21)) and then drives back to the depot  (pull-in arc (s21, k5)) to rest (waiting arc (k5, k6) and(k6, k7) ) and then drives out from depot  (pull-out arc (k7, s43)) to station s4 to perform another timetabled trip (task arc (s43, s33)) and finally drives back to depot  (pull-in arc (s33, k12)).There are two pull-out arcs and two pull-in arcs in the circuit, but the starting pull-out arc and the ending pull-in arc of the circuit are unique.The goal is to find a series of desirable circuits with minimum total cost.The number of circuits that start from the depot vertex is equal to the number of buses used in the scheduling system.

Network Modeling.
The key issue of the time-space network is to reduce the number of deadhead arcs.Most deadhead arcs can be removed as each of them can be represented by another deadhead arc and a waiting arc without reducing the number of compatible task pairs.For a large-scale system, the number of compatible task pairs is tremendous.Therefore, it is not desirable to connect all compatible task pairs by deadhead arcs beforehand and then reduce them.Instead, the two steps are carried out simultaneously.
For a depot  ∈ , the network is constructed with four steps.
Step 1 (generation of terminal station vertices and task arcs).Generate task arcs and corresponding terminal station vertices according to the known bus timetable.Every task arc has a start terminal station vertex and an end terminal station vertex.Every terminal station has several terminal station vertices that are arranged in a chronological order.
Step 2 (generation of depot vertices, pull-in arcs, and pull-out arcs).A time line connected by depot vertices that represents all possible arrival and departure events at the depot is built.
The pull-out arcs are constructed by connecting the departure depot vertices to the corresponding start terminal station vertices (generated in Step 1).Similarly, the pull-in arcs are constructed by connecting the relative end terminal station vertices (generated in Step 1) to the arrival depot vertices.
Step 3 (generation of waiting arcs).For a given terminal station, the waiting arcs outside the depot are constructed by connecting all of its terminal station vertices.For the depot, the waiting arcs inside the depot are constructed by connecting all of its depot vertices.
Step 4 (generation and reduction of deadhead arcs).Only those deadhead arcs that correlate to the first matches and the nearest matches are introduced into the model.Moreover, a very long deadhead arc (i,j), starting from terminal station vertex i at time t i and ending at terminal station vertex j at time t j , should also be removed if t j -t i >t, where t is the duration time of pull-in arc (from i to the depot) plus the duration time of pull-out arc (from the depot to j ).These techniques result in a dramatic reduction in the number of deadhead arcs.

Model Formulation
The MDVSP in this study is formulated as a set-partitioning model calling for the determination of a collection of circuits with minimum cost, which serves each task once and satisfies strict departure-duration restrictions.
Let K be the set of all depots.v k denotes the number of buses in depot k,  ∈ .Λ k is the set of all task arcs in the time-space network of depot k.Ω k represents the set of all possible circuits, each corresponding to a possible bus tour, which starts with a pull-out arc and ends up with a pull-in arc.Every circuit  ∈ Ω  has an associated cost c p .Apart from the variable operating cost (cost of different arcs covered by p), c p also includes fixed cost for the required bus.For every task arc i∈Λ k , let a ip be a binary coefficient that takes value 1 if task arc i is covered by circuit p and takes value 0 otherwise.For every circuit p∈Ω k , the binary decision variable x p takes value 1 if and only if the circuit p is selected in the optimal solution.Otherwise, x p takes value 0. The circuit p consists of a feasible sequence of different kinds of arcs chained with each other.The pull-out and pull-in arcs appear in pairs.Γ  is the set of pairs of pull-out and pull-in arcs.For a pair of pullout and pull-in arcs j, j∈ Γ  ,    denotes the start time of the pull-out arc and    denotes the end time of the pull-in arc.t 0 represents the working time limits of bus crews.The model is as follows. Minimize subject to: Mathematical Problems in Engineering Objective functions (7) serve to minimize the total scheduling cost of a bus company.Constraint (8) ensures that every timetabled trip is covered by exactly one bus.Constraint (9) imposes that the number of buses selected to run the timetabled trips should be no more than the total number in depot k.Constraint (10) imposes that the departure-duration of every selected bus should be no longer than the working time limits of bus crews.Constraint (11) defines the binary decision variables to depict whether a bus tour p is selected or not.

Solution Algorithm
In this section, an improved LNS algorithm based on the branch-and-price algorithm is constructed to solve the largescale MDVSP.We first present the framework of the algorithm before introducing the details of solving the restricted master problem and subproblem.

The Improved Large Neighborhood Search Algorithm.
The MDVSP explored in this study is an extremely large-scale Integer Programming (IP) problem.For one depot k, the number of binary decision variables is equal to the number of all possible circuits in Ω k .As the size of the circuit set (Ω k ,  ∈ ) is very large, we do not compute the set of all possible circuits when solving the IP problem.Traditional exact algorithms are not applicable.On the framework of LNS, an improved heuristic algorithm that is both accurate and efficient is put forward.As there are departure-duration restrictions in the explored large-scale MDVSP, the SPFA [26] is modified by adding an embedded preliminary exploring tactic to solve the subproblem.
The improved LNS algorithm adheres to the following steps.
Step 1 (obtain initial feasible solution).Generate an initial feasible solution S 0 by insertion heuristic.The reader is referred to [27] for detailed analysis of the insertion heuristic.
Step 2 (initialize parameters).Let parameters iter=0 and S best =S 0 , where iter denotes the number of iterations and S best denotes the best solution set at present.r represents the number of bus circuits selected from the best solution set S best at present.maxIter is the maximum number of iterations.
Step 3 (generate neighborhood).Generate a neighborhood by selecting r bus circuits from the best solution set at present S best .Select r bus circuits from the set S best .The set of r bus circuits is termed as S r .Generate a new problem using all tasks (timetabled trips) contained in the r bus circuits.This new problem has a smaller scale which will take much less time to be solved.The solution space of this new problem is the neighborhood.As to the selecting strategy of the r bus circuits, randomly select r bus circuits from those circuits that were selected less often in set S best .
Step 4 (search neighborhood).Search the neighborhood generated in Step 3 by branch-and-price and we get a better feasible solution of the new problem, termed as    .Reinsert    into set S best instead of S r .We obtain a feasible solution   of the original problem.
Step 7 (the end).S best is the best solution that can be found.
Figure 3 shows the flowchart of the improved LNS algorithm.The main idea of this proposed algorithm is to find an initial feasible solution S 0 at first and then improve this feasible solution through iteration until certain iteration stopping conditions are satisfied and we obtain the optimal solution.The initial feasible solution and the optimal solution are both a subset of (Ω k ,  ∈ ).As to how to improve this feasible solution S 0 , the main idea is to optimize a part of S 0 (r bus circuits in S 0 ) to get a better solution and repeat this operation again and again until certain iteration stopping conditions are satisfied.This operation is the process of generating a neighborhood and then searching the neighborhood.
The branch-and-price algorithm is based on the framework of the branch-and-bound algorithm.As the linear programming relaxation of the set-partitioning model is typically very tight [13], the original integer programming problem is relaxed into a linear programming (LP) problem (master problem) at each node in the branch-and-bound tree.
For a general LP problem (master problem): subject to: Q is the set of all columns and q is a column in Q.The simplex method can be used to solve this problem if the number of columns in Q is not large.Let  be the corresponding dual variable.In each iteration of the simplex method, if there's a column q whose   −   < 0, we should select   :   =  ∈ (  −   ).However, generally the number of columns in Q is very large, so we can solve this LP relaxation (master problem) using the column generation approach by first solving the LP relaxation over a subset of the possible columns.We refer to a master problem with only a subset of the possible columns as a restricted master problem.Additional columns can be generated as needed for the restricted master problem by solving the subproblem. (15) subject to: Let  * be the solution of this restricted master problem and  * be the corresponding dual variable.The subproblem can be expressed as  * =  ∈ (  −  *   ).For the MDVSP in this paper, the subproblem is the shortest path problem with departure-duration restrictions.Sections 4.2 and 4.3 elaborate the details of branch-and-price techniques for neighborhood searching.

Restricted Master Problem.
At each node in the branchand-bound tree, the original integer programming problem is relaxed into a linear programming (LP) problem (master problem).Relax the integer variable x p into a continuous variable between 0 and 1.Let Ω k 0 be the subset of Ω k .When applying column generation to solve the relaxed version of the original model in this paper, the master problem is divided into two parts, a restricted master problem and a subproblem.
ILOG CPLEX or Lingo can be applied to solve the DRMP above and obtain the value of decision variables ( i and  k ).The solution of the dual problem ( i and  k ) in this paper is obtained by applying ILOG CPLEX to the DRMP.Update the cost of arcs in the network according to the value of  i and  k to solve the subproblem.

Subproblem.
For the MDVSP in this paper, the subproblem is the shortest path problem with departure-duration restrictions.As there is a possibility of negative weighted arcs, the SPFA [26] is applied to solve the subproblem.In essence, SPFA is a Bellman-Ford algorithm with queue optimization.
As there are departure-duration restrictions in the MDVSP, the SPFA is modified by embedding a preliminary exploring tactic.Specifically, check if the vertex meets the departure-duration restriction before adding it into the queue.If the restriction is satisfied, add the vertex into the queue.Otherwise, do not add it.The depot vertices do not need to be checked by the preliminary exploring tactic.Only terminal station vertices have to be checked.Let u be the terminal station vertex that has to be checked by the preliminary exploring tactic.[] denotes the departureduration time of a bus tour that starts from a depot vertex and ends at vertex u. [] represents the corresponding tour cost.[] is the vertex just before u on the bus tour with minimal tour cost among all bus tours that start from a depot vertex and end at vertex u.For a terminal station vertex u, traverse all arcs (including task arcs and waiting arcs) starting from u.If it is a task arc and []+t+td≤  0 , where  is the trip time of the task arc,  is the time from the end of the task arc to a depot, and  0 is the known working time limits, continue the SPFA.Otherwise, vertex u cannot be added to the queue.If it is a waiting arc, continue searching until a task arc appears.
The flowchart of the modified SPFA is shown in Figure 4.

Case Study
In this section, a real-life case in China with 3 depots, 8 bus lines, and 930 timetabled trips is selected to evaluate the performance of the proposed methodology.Some basic information about the real-life case is shown in the "Supplementary Materials" (available here).There are 96 buses in total for the real-life case with 52 buses in Depot 1, 23 buses in Depot 2, and 21 buses in Depot 3. The proposed algorithm is implemented by c# language in Visual Studio 2005.ILOG CPLEX is used to solve the DRMP.In this paper, let r=20.The program runs on a computer with CPU Intel T7500 2.20GHz 2G.The optimal schedule scheme for the real-life case is shown as Table S5 in the "Supplementary Materials".87 buses are used to perform the 930 timetabled trips with 48 buses in Depot 1 performing 511 timetabled trips, 21 buses in Depot 2 performing 207 timetabled trips, and 18 buses in Depot 3 performing 212 timetabled trips, respectively.
To highlight the efficiency of the improved LNS algorithm, the computational result is compared with the result using the branch-and-price algorithm.The branch-and-price algorithm is also implemented by c# language in Visual Studio 2005 and runs on a computer with CPU Intel T7500 2.20GHz 2G.The comparison of the two results is shown in Figure 5.
As shown in Figure 5, both of the algorithms have a fast convergence speed at first and then slow down.However, compared to the improved LNS algorithm, the number of feasible solutions from the branch-and-price algorithm is much smaller.Moreover, it takes much longer to obtain the first feasible solution for branch-and-price.The total cost of the improved LNS algorithm is 299,672 which is slightly more than that of branch-and-price (299,017).This result does make sense as the LNS algorithm optimizes the problem locally while the branch-and-price optimizes the problem globally.The difference in solution quality is due to the difference between local optimization and global optimization.Therefore, the branch-and-price obtains higher quality solutions in the real-life case study.But, as we can see from the results, the difference in the total cost is very small, which can be ignored.However, the improved LNS algorithm begins to converge at 34 minutes, which is much shorter than 104 minutes using branch-and-price.In summary, the improved LNS algorithm can ensure not only the quality of the solution but above all the computational efficiency.On that account, the proposed LNS algorithm is of great significance for large-scale MDVSP, which attaches great importance to the algorithmic efficiency.
In Figure 5 the total cost of the improved LNS algorithm has two hops at 12 minutes and 31 minutes during the iteration process.The reason probably lies in the fixed cost of buses.Compared to other costs like   ,   in Section 2, the fixed cost of a bus is relatively larger.The objective function's value may change significantly as the number of buses needed decreases during the optimization process.

Sensitivity Analysis
Although the proposed methodology is tested using a real-life case, extra computational experiments should be designed to further validate whether the algorithm is appropriate for large-scale systems.The most effective method for this purpose is to compare the performance of the improved LNS algorithm and branch-and-price in terms of quality and efficiency, using test instances with different number of timetabled trips.Test instances in this section are generated in a 6km * 6km square using the same methods as Carpaneto et al. [28].Test instances with n=500, 1000, 1500, and 2000 timetabled trips are generated.For each n, 10 test instances have been generated as done by Carpaneto et al. [28].To be comparable, let the number of depots K=3 (which is the same as the real-life case in Section 5).Quality Gap (QG) is used to compare the solution quality of the two algorithms and Efficiency Ratio (ER) is used to compare the computational efficiency of the two algorithms.QG and ER are defined as follows: where   and   are the total cost and computational time of the best solution obtained by branch and-price.  and   are the total cost and computational time of the best solution obtained by the improved LNS algorithm.
All the test instances are run on a computer with CPU Intel T7500 2.20GHz 2G.The computational results of the test instances with different number of timetabled trips are presented in Table 1.As seen in Table 1, the total cost of the improved LNS algorithm is slightly more than that of branch-and-price algorithm for each test instance of different trip size (500, 1000, 1500, and 2000).But, as we can see from the Quality Gap, the difference in the total cost is very small, which can be ignored.However, the computational time of the improved LNS algorithm is much less than that of branch-and-price algorithm for each test instance of different trip size (500, 1000, 1500, and 2000) which indicates that the improved LNS algorithm is much more efficient than branch-and-price, and  the advantages of algorithmic efficiency continue to expand as the number of timetabled trips increases.
To test whether there are statistically significant differences between the improved LNS algorithm and the branchand-price algorithm in total cost and computational time, "Wilcoxon Signed Ranks Test" in SPSS is applied to the data in Table 1.The results are shown in Table 2.The variable name in Table 2 consists of two parts.Take "BP Cost", for example, "BP" represents the method name (branch-and-price) and "Cost" represents the total cost.
As we can see from the Ranks table in Table 2, LNS algorithm had a higher total cost but a lower computational time than the branch-and-price algorithm for all test instances.By examining the final Test Statistics table in Table 2, we can discover whether there are statistically significant differences between the improved LNS algorithm and the branch-and-price algorithm in total cost and computational time.We are looking for the "Asymp.Sig.(2-tailed)" value, which in this case is 0.000, illustrating that the observed differences are statistically significant.However, compared with the significant improvement in algorithmic efficiency brought by the LNS algorithm, the increase in total cost is very small, which can be ignored.
In order to better show the results, the minimum value, lower-quartile value, median value, upper-quartile value, maximum value, and average value of Quality Gap and Efficiency Ratio over 10 test instances for each n (n=500, 1000, 1500, 2000) are calculated, as shown in Tables 3 and 4.
As seen in Tables 3 and 4, there is not a big change in Quality Gap or Efficiency Ratio for different test instances of each n (n=500, 1000, 1500, and 2000).The Quality Gap is approximately 0.35%, 0.38%, 0.63%, and 0.93%, respectively, Based on Tables 1 and 2 and the computational results of Quality Gap and Efficiency Ratio presented in Tables 3 and  4, the findings suggest that the improved LNS algorithm is effective in obtaining faster solutions without deteriorating the quality, especially for large-scale instances.

Conclusions and Future Work
To solve this large-scale MDVSP, the only feasible approach is to downsize the problem and improve the algorithm simultaneously.Crews working overtime increase safety risks for themselves and others, so buses should return to the depot and change shifts when the crews reach their working time limits.Therefore, an improved set-partitioning model with departure-duration restrictions was established based on the time-space network, calling for the determination of a collection of circuits with minimum cost.In order to solve the MDVSP model, an improved LNS algorithm was proposed.To deal with departure-duration restrictions, the SPFA was modified by embedding a preliminary exploring tactic to solve the subproblem.A case study in China was presented to evaluate the performance of the proposed algorithm.The results indicate that the improved LNS algorithm, which encompasses attributes of both fast computing speed and high quality solutions, is especially appropriate for an extremely large-scale MDVSP.The sensitivity analysis in Section 6 further verifies the effectiveness of the proposed methodology to solve a large-scale MDVSP.
In this study, only bus scheduling is considered, which is just one part of the whole bus operating planning process.Future study should be directed at researching all four parts (i.e., bus network design, timetable design, bus scheduling, and crew scheduling) and their integration to continue developing a more stable and sophisticated bus scheduling plan.In addition, in the proposed LNS algorithm, the selecting strategy of the r bus circuits and the value of parameter r affect the computational efficiency of the proposed algorithm.A good selecting strategy and r value will further improve the LNS algorithm.The selecting strategy is a recommended future research topic.How to scientifically select the value of parameter r also needs further research.closing time, and the number of timetabled trips of each bus line, (3) travel time between the depots and the origin (destination) stations of each bus line, (4) travel time between different origin and destination stations of each bus line, (5) timetable information of each bus line, and (6) the optimal schedule scheme for the real-life case are shown as Figure S1, Table S1, Table S2, Table S3, Table S4, and Table S5 in the "Supplementary Materials".In Table S2 and Table S3, " i-O" represents the origin station of Line i (upward) and " i-D" represents the destination station of Line i (upward).The origin station of an upward bus line is the destination station of its corresponding downward bus line and the destination station of an upward bus line is the origin station of its corresponding downward bus line.Some lines have the same origin station or destination station; for example, Line 1 (upward), Line 2 (upward), and Line 6 (upward) have the same destination station.(Supplementary Materials)
solution S 0 S ＜？ＭＮ =S  S  <S ＜？ＭＮ ?Select r bus circuits from S ＜？ＭＮ to generate Let iter=0, S ＜？ＭＮ =S 0 neighbourhood and get a feasible solution S  Set the value of r and
Time < LNS Time b. BP Time > LNS Time c. BP Time = LNS Time Test Statistics a BP Time-LNS Time Z -5.511 b Asymp.Sig.(2-tailed) .000a. Wilcoxon Signed Ranks Test b.Based on negative ranks.

Table 1 :
Computational results for test instances with different number of timetabled trips.

Table 2 :
Results of Wilcoxon Signed Ranks Test.

Table 4 :
Efficiency ratio (ER).test instances with different number of timetabled trips (500, 1000, 1500, and 2000), while the Efficiency Ratio reaches up to 2.89, 2.98, 3.65, and 3.79, respectively, for test instances with different number of timetabled trips (500, 1000, 1500, and 2000).The small Quality Gap demonstrates the effectiveness of the improved LNS algorithm to produce high quality solutions.Meanwhile, the large Efficiency Ratio clearly highlights the excellent computational efficiency of the proposed heuristic. for