Exact and Metaheuristic Approaches for a Bi-Objective School Bus Scheduling Problem

As a class of hard combinatorial optimization problems, the school bus routing problem has received considerable attention in the last decades. For a multi-school system, given the bus trips for each school, the school bus scheduling problem aims at optimizing bus schedules to serve all the trips within the school time windows. In this paper, we propose two approaches for solving the bi-objective school bus scheduling problem: an exact method of mixed integer programming (MIP) and a metaheuristic method which combines simulated annealing with local search. We develop MIP formulations for homogenous and heterogeneous fleet problems respectively and solve the models by MIP solver CPLEX. The bus type-based formulation for heterogeneous fleet problem reduces the model complexity in terms of the number of decision variables and constraints. The metaheuristic method is a two-stage framework for minimizing the number of buses to be used as well as the total travel distance of buses. We evaluate the proposed MIP and the metaheuristic method on two benchmark datasets, showing that on both instances, our metaheuristic method significantly outperforms the respective state-of-the-art methods.


Introduction
Due to the large scale of school merge and closures in China since the late 1990s, more and more primary and secondary schools have begun providing school bus service. According to the authors' survey, in 2013, there are already thousands of school buses in service in cities such as Shanghai, Guangzhou, Shenzhen, Changchun, Xiamen and Qingdao. However, most of the school buses are operated by schools respectively. Recognizing the importance of school bus safety and operation efficiency, some local authorities, such as Binzhou, Chifeng, Wuli, Youxian and Zhuzhou, have proposed to provide regional level bus service for all schools. The Ministry of Education of China has launched a statewide school bus pilot program in six counties in 2011, aiming to explore the best practice of school bus operations in China. For local government and transportation company, school bus route planning is already a difficult task; still, it is crucial to serve with quality, equality and efficiency.
problem respectively. Section 4 shows the computation results for benchmark instances. The optimization performance of exact method and metaheuristic algorithm are compared in detail. The last section offers the research findings and conclusions.

Problem definition
In SBRP, a number of prerequisites, such as a set of schools, a set of student stops, a set of depots, a fleet of school buses and their geographic background, must be available for school bus route planning. The SBRP dataset should include both geographic and attribute information. The bus travel distances and times between schools, stops and depots can be prepared by geospatial analysis of the transportation network in GIS. The objective of school bus route planning is to minimize the number of school buses and the total bus travel distance while satisfying service qualities such as the student maximum riding time in a bus, the earliest and latest arrival times to their schools, and the maximum bus capacity.
Fig 1 illustrates two bus route plans for an example of school district. There are 2 schools and 8 student stops in the school district (Fig 1a). The student stops are filled with the same color as their schools. For school bus service in this district, Fig 1b shows a single-school bus route plan, in which there are two trips for school s 1 and s 2 respectively. Each trip consists of the depot, a sequence of bus stops and their destination school. One bus starts from the depot, picks up students of each bus stop in its trip, delivers them to their school, and then returns to the depot. When two trips can be served by the same bus in sequence without violating service constraints, they are assigned to the same bus (Fig 1c). This is a school bus scheduling problem, which may reduce the number of buses required. In Fig 1c, the dotted line represents the empty travel trip between two connected trips.
This paper focuses on the bus scheduling problem. In following sections, the MIP models for school bus scheduling problems are formulated in detail. Note that the bus trip in this paper refers to a bus route for a single school. The bus scheduling is to assign buses to trips, thus one bus can serve one or more bus trips.

MIP formulation for homogeneous fleet school bus scheduling
If the bus trips for each school are given, the bus scheduling problem is to assign buses to trips, and find the optimal solution with minimum number of buses and minimum bus travel distance. Let V' be the set of single-school bus trips, and the number of trips is n. Each rip is composed of a sequence of bus stops and their designation school. The travel distance,d ij , from trip i to trip j, can be determined by the travel time from the designated school of trip i, to the first  Table 1 shows the parameters and decision variables for bus scheduling with a homogeneous fleet of buses.
The bi-objective MIP formulation for homogeneous fleet school bus scheduling problem (MIPHO-B) is shown as follows: Subject to : eðsðiÞÞ y i lðsðiÞÞ; 8i 2 V 0 ð4Þ y i þ t ij þ p j À y j ð1 À x ij ÞM; 8i 2 V; j 2 V 0 ji 6 ¼ j ð5Þ x ij 2 f0; 1g; 8i; j 2 V 0 ji 6 ¼ j ð6Þ The objective Eq (1) is to minimize the number of buses and total travel distance. The first part of Eq (1) denotes the number of buses. The parameter M 0 is a big number that ensures the number of buses is as the primary objective and minimizes the two objectives in a lexicographic manner. Eqs (2) and (3) ensure that each trip is served by exactly one bus, and the bus must leave to another trip or the depot. Eq (4) assure that the bus serving trip i should arrive to school s(i) at any time within the school time window. Eq (5) state the accumulation of bus travel time:(y i + p j + t ij −y j )x ij = 0. If the arc (i,j) is traversed by a bus (x ij = 1), then y i +p j +t ij = y j . The non-linear constraints can be transformed to linear inequalities Eq (5) by using the big M parameter. Eqs (6) and (7) are the constraints on decision variables.
The travel time from i to j, i, j 2 V |i 6 ¼ j e 0 ,l 0 The earliest and latest departure time for buses from the depot Decision variables x ij If a bus serves trip i and leaves to serve trip j immediately, x ij = 1; otherwise x ij = 0, i, j 2 V |i 6 ¼ j

MIP formulation for heterogeneous fleet school bus scheduling
Unlike the assumption of a set of identical buses in the homogenous fleet model, a fleet of buses usually have different capacities and purchase cost in the heterogeneous fleet problem. There are two approaches of defining the three-index decision variables x ijk in problem formulation, where k presents a vehicle [17,28] or a type of vehicles [29]. Kim et al. developed a busbased MIP model (MIPHE) for heterogeneous fleet school bus scheduling problem [25], and solved a set of problem instances by CPLEX. We rewrite the model by changing the bus-based decision variables to the bus type-based variables, aiming to reduce the model complexity. In the bus-based formulation, k in the decision variables x ijk denotes the bus k to serve arc (i,j) in which i and j denotes the single bus trip. But, in the bus type-based formulation, the k denotes the bus of type k to serve arc (i,j). In addition, the fixed cost for each type of buses is also introduced in our objective function, which is different from that defined in MIPHE. Most parameters of the formulation are the same as the homogeneous problem defined in Table 1; additional parameters and decision variables are shown in Table 2.
The bi-objective MIP formulation for heterogeneous fleet school bus scheduling problem (MIPHE-B) is shown as follows: Minimize : Subject to : x ijk 2 f0; 1g; 8i; j 2 Vji 6 ¼ j; k 2 K ð16Þ y ik 2 f0; 1; 2; . . .g; The objective Eq (8) is to minimize the total fixed cost of buses and total travel distance. Eqs (9) and (10) ensure that each trip is served by only one bus, and the bus shall then leave to another trip or the depot. Eq (11) ensure that the number of used buses of type k does not exceed the total number of type k buses. Eqs (12) and (13) assure that a bus of type k must have enough capacity to serve a trip. Eq (14) assure that the bus of type k serving route i should arrive to school s(i) at any time within the school time window. Eq (15) state the accumulation of bus travel time: (y ik +p j +t ij -y ik )x ijk = 0. Eqs (16) and (17) are the constraints on decision variables. The bus type-based formulation for the school bus scheduling problem can reduce the model complexity significantly. Comparing to the bus based model MIPHE, the number of decision variables is reduced from (n 2 +n)m to (n 2 +n)k, where n, m and k represent the number of bus stops, the number of buses and the number of bus types. For an example of problem instance with 50 bus stops and a fleet of 20 buses, we need to define 51,000 decision variables in the bus-based formulation. If the buses are categorized into 3 types, we need to define 7,650 decision variables in bus type-based formulation, where 85% of decision variables are reduced. The parameters F k in the objective function encourage the model solver to select small buses with smaller fixed costs. In model building for a problem instance, the model can be simplified by introducing parameters E ik . If E ik E jk = 0 (or E ik = 0), the decision variable x ijk = 0 (or y ik = 0) can be deleted from the model, and the related elements in the objective function and constraints can also be deleted. In fact, all the Eqs (12) and (13) can be removed from the model. In addition, if e(s(i))+t ij +p j >l(s(j)), the connection between trips i and j is infeasible. In such a case, the decision variable x ijk = 0; thus the related elements in the objective function and constraints can also be deleted, and some Eq (15) will also be deleted.
Since the difference of capacity and cost between a fleet of buses is considered in the heterogeneous fleet problem, it is more complex than the homogenous fleet problem. Technically, the MIP model for the homogeneous fleet problem is a special case of the heterogeneous fleet problem, where all the buses are assumed to be identical. In addition. Compared with the bus-based model MIPHE, the bus type-based model MIPHE-B reduces the number of decision variables significantly; thus, its computational complexity can be reduced.

Algorithm framework
In this paper, the multi-school SBRP is decomposed into single school bus routing problem and school bus scheduling problem. The former problem is equivalent to open CVRP with a constraint of maximum riding time of students. The latter can be converted to VRPTW by treating each trip as a virtual bus stop [25]. Based on the neighborhood search operators for VRP, we propose a two-stage algorithm framework for both problems in the following steps: Step 1 reads the data of a SBRP instance.
Step 2 generates a set of feasible routes as the initial solution. The feasible routes must meet specific constraints.
Step 3 is to minimize the number of routes using neighborhood operators coupled with simulated annealing. Since the multi-objective function in step 3 encourages reducing the number of bus routes, additional step 4 is used to minimize the total travel distance. The output solution in step 3 is fed as the initial solution of step 4.
There are two general objectives for SBRP in this research: minimizing the number of buses and minimizing the total travel distance. Since the capital cost per bus is significantly higher than the travel cost per distance unit [30], the number of routes is of higher priority than the travel distance. However, it is very difficult to obtain a solution that all the objectives are optimal simultaneously. The objective function of total travel distance commonly used in VRP may lead the search to the solutions with a small travel distance and makes it impossible to remove routes [31]. Therefore, we use a lexicographic multi-objective function to evaluate the neighbourhood solutions in step 3: In the function, the |R|, Dist(R)and |r| 2 are the number of routes, the total length of all routes, and the sum of squared route stop numbers respectively. Given the parameters α ) β ) γ, the neighbourhood moves prefer to delete some routes if possible. In step 4, only the function of Dist(R) is used to minimize the total travel distance. After the construction of an initial route solution by a simple heuristic algorithm, the key steps of this algorithm are to improve the solution by exploring the neighborhood solutions. Local search is a metaheuristic method for solving combinatorial optimization. Various local search algorithms iteratively improves the solution of a problem by applying local changes to current solution until optimal solution is deemed to be found or a stop criterion is met. In the process of searching a better solution, once a feasible neighborhood solution is identified and its objective is better than the current solution, the current solution will be replaced by the new one. Local search has proven to be an effective technique for generating good solutions to instances of the classical VRP as well as variants such as VRPTW [32].
A general local search operator exchanges μ nodes from one route with λ nodes from another. The local search operators such as one-point move, two-point move and 2-opt move can work on a single route or different routes, while cross-exchange is an inter-route operators that only work on two routes simultaneously [33]. One-point move attempts to relocate an existing node to a new position, which can reduce the number of routes and create more opportunities for other operators. Two-point move swaps two nodes to reduce route distance. 2-opt move tries to improve a solution by removing two edges from the solution and replace them with two new edges. Cross-exchange move tries to remove two edges respectively from two different routes and replace them with four new edges. The repeat of the four local search operators will improve the route solution gradually.
The local search algorithm has drawbacks: the quality of solution is closely related with the initial solution and the spatiotemporal structure of the problem instance. Meanwhile, the procedure is very easily trapped into local optima. In general, such a local optimum is not a global optimum [34]. Even worse, there is usually no guarantee that the value of the objective function at an arbitrary local optimum is close to the optimal value. Some methods have been proposed to solve the problem of local optimum, such as simulated annealing and tabu search. In this paper, the technique of simulated annealing [35] is combined with local search operators. In the local search operations, we accept some solutions that make route length longer by simulated annealing, but may escape from local optima and realize the search in depth and diversification. The general local search procedure coupled with simulated annealing is designed in the following steps: This procedure can be implemented for school bus routing and bus scheduling for the purpose of minimizing the number of routes and minimizing the total travel distance. Line 4 permutes the nodes randomly, which will change the node searching sequence in each iteration for solution diversification. From line 5 to line 8 in the loop, one-point move search, two-point move search, 2-opt move search, and cross-exchange search shall be executed sequentially. In each local search function, according to feasibility check, objective evaluation and acceptance criterion; the neighborhood solution will either be accepted as the current best solution so far or rejected. In line (9) the temperature drops, and the possibility of accepting a new solution with longer travel distance in next loops will be decreased gradually.

Local search operators
The route solution sol is gradually improved by four local search operators: one-point move, two-point move, 2-opt move and cross-exchange move. Algorithm 2 shows the details of the one-point move search. The search procedures for other three operators are similar to onepoint move search. Input: the SBRP instance (S), the sequential or randomized node list (perm), and the criteria to accept a new solution (rule). Note that the instance S includes the current best route solution (sol best ), the current temperature (t) and other information

}
Output: the current best solution (sol best ) For each node p, the algorithm constructs its neighborhood list (Nlist) at first (line 2). For each neighbor node p Ã in the neighborhood list, if the node p can be inserted after node p Ã without violation of the constraints, the objective function for this feasible move will be evaluated. Otherwise, the operator will try another neighbor position (line 4). If the neighborhood solution (line 5) is better than the current best solution sol best or, at least acceptable with the simulated annealing criteria (line 6), the best solution will be updated according to acceptance criterion (rule): first acceptance (FA) or best acceptance (BA) (line [7][8][9][10][11][12]. The function is_better(sol best , sol Ã ) evaluates the solution quality of sol Ã compared with that current best solution sol best ; The function is_better2(sol Ã best , sol Ã ) is used to select the best neighborhood solution. For minimizing the number of routes, the three elements in objective Eq (18) are in a lexicographic order.

Solvers for school bus routing and scheduling
Based on the algorithm framework and local search operators discussed in Section 3.1 and 3.2, we developed three solvers for single-school SBRP, homogeneous fleet school bus scheduling problem and heterogeneous fleet school bus scheduling problem respectively. The solvers are used to plan bus routes. Firstly, the bus trips for each school in a multi-school system is generated by the single-school solver. Secondly, the buses are assigned to the trips by using the homogeneous fleet or heterogeneous fleet school bus scheduling solver. The output data from the first solver is fed into the second solver. In each solver, a simple savings-based [36] heuristic algorithm is used to construct an initial solution. A two-stage metaheuristic algorithm is then used to improve the solution. The first stage minimizes the number of routes, which reduces the total number of buses needed for students. The second stage optimizes the total bus travel distance, which will reduce the total bus travel cost. The procedures for the two stages are designed in the same manner, except for the evaluation of solution objectives.
The VRP with a heterogeneous fleet is considered much harder than corresponding problems with a homogeneous fleet [37]. The heterogeneous fleet school bus scheduling problem is also much more difficult to solve than the homogeneous fleet problem. For initial solution construction, we divide the heterogeneous fleet problem into multiple homogeneous fleet subproblems according to the types of buses. First, the trips and the types of buses are sorted in an ascending order. Second, the trips are divided into multiple groups according to the capacities of buses; a type of buses with minimum capacity is assigned to each group of trips as well. Third, the saving-based heuristic algorithm is used to construct the route solution for each group of trips. In case the number of one type of buses is not enough to serve the sub-problem, the types of buses for all routes need to be adjusted appropriately. Fourth, the route solutions of each group are merged into an initial route solution for the whole heterogeneous fleet problem. The initial solution is then improved by the two-stage algorithm. Comparing with the neighborhood point moves in solving homogeneous fleet problem, it is needed to check the bus types for any possible point moves.
The solver programs are implemented by C++ programming and compiled as 32-bit command-line executables with parameter setting options. The parameters includes the maximum number of loops (maxloop), the initial temperature for simulated annealing(t 0 ), the cooling rate of temperature(α), the length of neighborhood list (NList), the acceptance criteria for solution updating (rule), and the randomized node search strategy (random).

Results for HOMO and HETERO instances
The HOMO and HETERO benchmark instances are collected from literature [25]. There are 15 instances for homogeneous fleet problem and heterogeneous fleet problem respectively. Table 3 shows the size of all instances. For each instance, the school data includes the school ID, the spatial coordinates x and y and the school time window. The bus trips for each school are given, which include the trip ID, the coordinates of the first stop in the trip, destination school ID, and the bus service time. The bus data for heterogeneous fleet instances include the bus ID and the bus capacity.
We build models for all the instances according to the MIP formulations proposed in Section 2.2 and 2.3. The model parameters are the same as literature [25]: the travel distance d ij is an integer number calculated by Euclidean distance with a truncated fraction; the bus speed is assumed to be 1 distance unit per time unit, and thus t ij = d ij . In the objective function, the travel distance from the bus depot to the first trip and from the last trip to the depot are not considered in the total distance travelled by all the buses. In the model MIPHE-B, the bus type data are counted from the bus data. There are three types of buses in the heterogeneous fleet instances; and the fixed costs for each type of buses are assumed to be 100000, 90000 and 80000 for buses with capacity of 70, 50 and 40 respectively. The parameters M and M 0 are set to 10 9 and 10 5 respectively. All the models are solved by the MIP solver IBM ILOG CPLEX 12.6 (64bits). CPLEX is installed on a personal computer with the following hardware configuration: Intel(R) Core2 Duo E7600 3.06G CPU, 4GB RAM and 64-bit Windows 7 operating system. The parameters for CPLEX solver are set to its default values, except that the maximum computation time is set to 3600s and the MIPGap is set to 10 −12 . Tables 4 and 5 show the CPLEX results of the homogeneous fleet and heterogeneous fleet instances respectively. The columns N, D and T provide the number of buses required, the deadhead bus travel distance, and the CPU time that are used by the CPLEX solver. The deadhead bus travel distance is defined as the idle distance covered by the bus between the destination school of the first trip to the first bus stop of the next trip without carrying any students, when a bus can serve the two trips in sequence. The column Gap indicates the CPLEX solution status: optimal or the optimality gap. The instances of homo1~homo60 can be solved with optimal solutions; while the instances of homo70~homo100 can be solved with suboptimal solutions. Since the new version of CPLEX is used in our research, we get better solutions with the less number of buses for large instances of homo70~homo100. For the last two instances (homo90 and homo100), we get longer travel distances. However, we provide the solutions with less buses, which will reduce the total cost of school bus service. Note that our optimal solution for homo60 is different from that in literature [25]. Table 5 also shows the number of each type of buses used for the instances. The numbers in bracket denote the number of buses with capacity of 40, 50 and 70 respectively. The CPLEX results show that our model MIPHE-B for the heterogeneous fleet school bus scheduling problem outperforms the model MIPHE significantly. Since the number of buses is rather large compared to the types of buses, the bus-based formulation of MIPHE is computationally more complex than the bus type-based formulation. Thus, the CPLEX can only solve the small instances (hetero1~hetero10). Our bus type-based formulation MIPHE-B reduces the model complexity in terms of the number of decision variables, the number of constraints, and the model file size. Modeling the instances with MIPHE-B, CPLEX can solve larger instances (het-ero20~hetero50) optimally; the optimal problem size is increased from 58 (hetero10) trips to 286 trips (hetero50). For the large instances (hetero60~hetero100), the CPLEX can provide feasible solutions, and the results for some instances (hetero60~hetero80) are near to the optimal solutions. At the same time, since the difference between the fixed costs of buses is not considered in MIPHE model, the optimal solutions for a instance obtained from MIPHE and MIPHE-B may be different and cannot be compared directly. For example, the bus travel distances for hetero5 and hetero10 from MIPHE-B are longer than that from MIPHE. However, we provide the optimal configuration of the three types of buses. The parameter of fixed bus cost (F k ) in MIPHE-B model can guide the MIP solver to select smaller, lower cost buses as many as possible.
The heuristic results for HOMO and HETERO instances are listed in Tables 6 and 7 respectively. Since simulated annealing is a generic probabilistic metaheuristic, the solver proposed in Section 3 was executed 10 times over each instance. For both homogeneous fleet and heterogeneous fleet instances, the size of neighborhood list (NList) is set to n/2, where n is the number of bus trips. The initial temperature for simulated annealing is set to 5.0/n, and the cooling rate α is 0.995. The maximum number of loops (maxloop) is set to 200. The best, average and standard deviation of the number of buses and the total distance among the 10 solutions are listed in columns N best , N ave , N dev or D best , D ave , D dev , which are useful to verify the effectiveness of the algorithm. Table 6 shows the heuristic results for the homogeneous fleet instances. The columns N best , N ave and N dev represent the best, average and standard deviation of the number of buses among the 10 solutions respectively. The columns D best , D ave and D dev represent the best, average and standard deviation of the distance among the 10 solutions respectively, and the column N gap and D gap represents the percentage difference of the number of buses and distance between the best heuristic solution and the CPLEX solution for each instance in Table 4, respectively. The column T denotes the solver time in seconds. Compared with CPLEX, the heuristic approach to homogeneous fleet school bus scheduling can find the optimal solutions for small-size instances (homo1~homo10), sub-optimal solutions for medium-size instances (homo20~homo60), and better solutions for large-size instances (homo90 and homo100). In addition, compared with the solution results in literature [25], our SA-heuristic algorithm outperforms Kim's heuristic algorithm in terms of the number of buses and the total travel distance; for large instances homo60~homo100, on average, the number of buses and the total travel distance are reduced by 2.35% and 2.83% respectively. Table 7 shows the heuristic results for the heterogeneous fleet instances. Compared with the solutions in literature [25], we provide the number of buses by bus type. In addition, the objective differences between the heuristic solutions and CPLEX solutions are shown in column Gap1, Gap2 and Gap3, which denote the percentage gaps of the number of buses, the fixed cost of buses and the travel distance respectively. The solution comparison between exact and heuristic approach shown in Fig 2 is interesting. For small-size instances (hetero1~hetero5, the number of trips between 7 and 32), both the two approaches can provide optimal solutions in terms of the number of buses, the total fixed cost and the total travel distance. For medium-size instances (hetero10~hetero50, the number of trips between 58 and 286), CPLEX can find the optimal solutions, and our algorithm can find good enough results. For example, our algorithm can find the optimal number of buses on instances hetero10 and hetero20. For large-size instances (hetero60~hetero100, the number of trips between 342 and 562), CPLEX is able to find sub-optimal solution with small objective gaps for instance hetero60~hetero80, whereas high quality solutions for instances hetero90 and hetero100 are difficult to find. However, the Approaches for a Bi-Objective SBSP heuristic approach has obvious advantage in solving the larger instances in both the solution quality and computation time.

Results for RSRB and CSCB instances
Another set of sixteen benchmark instances are selected from literature [27]. The instances shown in Table 8 are classified into two groups: the random spatial distribution of schools and bus stops (RSRB) and the clustered distribution (CSCB). These instances have different numbers of students, bus stops and schools. The capacity of the vehicle is assumed to be 66, and its average speed is 20 miles per hour. For the purpose of comparison, the service time of picking up or delivering students at each stop is calculated by the regression model developed in [20], the same as in literature [27]. The maximum riding time (MRT) of students in a bus is set to 2700 seconds (45 minutes) or 5400 seconds (90 minutes). Since there is no bus depot information in the benchmark instances, we assume the bus depot is located at coordinates (105600, 105600). For school bus scheduling, the bus trips for each school are solved by the metaheuristic algorithm proposed in Section 3. The travel distance d ij between any two nodes is a float number in inches and is calculated by Manhattan distance. The travel time t ij between any two nodes is an integer number in seconds, it is estimated by the travel distance and the average bus speed, i.e. t ij d ij /29.333333. The pickup time at student stop is an integer number estimated by the formula 19+2.6 Ã n (n is the number of students at the stop); and the delivery service time at school stop is also an integer number, which is estimated by the formula 29+1.9 Ã n (n is the number of students served by the trip). The output trip data includes these features: the trip ID, the first student stop ID and its coordinates, the destination school ID and its coordinates and time windows, the bus travel distance, the trip service time, and the number of students served by the trip. The service time for each trip includes the bus travel time, the students' pickup time and unloading time at the destination school. The column #trips1 in Table 8 presents the total number of trips for each instance when the MRT is 2700 seconds; and column #trips2 presents the total number of trips for each instance when the MRT is 5400 seconds.
Based on the trips for each school, we build all the models for the instances according to the MIPHO-B formulations proposed in Section 2.2. The parameter M in objective function is set to 10 9 . Parameters such as bus travel distance and time are also calculated, same as the settings in single-school problem. The models are solved by CPLEX 12.6. The CPLEX parameters are set to default values, except that the maximum computation time is set to 3600s and the MIP-Gap is set to 10 −15 .
All the instances in Table 8 are also solved by our solver for homogeneous fleet school bus scheduling problem. The solver was executed 10 times over each instance. The initial temperature for simulated annealing is set to 1.2/n, where n is the number of bus trips of each instance. The cooling rate α is set to 0.99. The size of neighborhood list (NList) is set to n/2; and the maximum number of loops (maxloop) is 200.
Tables 9 and 10 exhibit the computational results for RSRB and CSCB instances respectively. The 2700 and 5400 in column Instance denote the students' maxmium riding time in seconds. The columns N, D and T denote the total number of buses, the total bus travel distance and the solver time in seconds. The column Gap indicates the solution status for CPLEX solver; the columns N gap and D gap refer to the percentage difference of the number of buses and the travel distance between CPLEX solutions and the heuristic solutions. For the metaheuristic solver, the columns N best , N dev , D best , D dev have the same meaning as the corresponding columns in Table 6.
There are several findings from Tables 9 and 10. First, the number of buses can be reduced significantly by bus scheduling. Comparing to the single-school based bus route plan, only about 30% of buses are needed for both RSRB and CSCB instances. Second, the homogeneous fleet bus scheduling problem, which is modeled as VRPTW without capacity constraints, can be efficiently solved by CPLEX with optimal solutions. Third, the local search algorithm coupling with simulated annealing is also effective to solve homogeneous fleet bus scheduling problem. The algorithm proposed in this research can also find the minimum number of buses for the benchmark instances except a few instances. Table 11 shows the comparison of our algorithm with the heuristic algorithm in [27]. Z Park and Z SA denote the number of required buses obtained from the two algorithms. The Gap column shows the percentage improvement of our algorithm, which is computed as (Z Park -Z SA )/ Z Park . On average, our route solutions reduce the number of required buses by 24.8% and 15.98% for RSRB and CSCB instances respectively. It is mainly due to Park's algorithm is a construction heuristic algorithm with simple improvement operators only.

Conclusions
For a multi-school system, the school bus scheduling problem seeks to optimize bus schedules to serve all the given trips considering the school time windows. Two approaches to solve the bi-objective school bus scheduling problem are discussed in this paper. This research manages to make two contributions: the rewriting of MIP formulation for heterogeneous-fleet school bus scheduling, and the design of a local search algorithm coupling with simulated annealing to minimize the number of buses and the total travel distance.
This research also reports three findings. First, the school bus scheduling problem can be efficiently solved with optimal or sub-optimal solution by MIP solver CPLEX. School bus scheduling problems is a class of combinatorial optimization problem, and is also NP-complete Approaches for a Bi-Objective SBSP problem. However, along with the great progress in integer programming, especially the introduction of hybrid heuristics in MIP solver design, it is possible to solve large-size real-world school bus scheduling problems by MIP optimizer. For homogenous fleet problem, the CSCB and RSRB instances with up to 745 trips and the HOMO instances with up to 372 trips are solved optimally. For heterogeneous fleet problem, the HETERO instances with up to 286 trips can be solved optimally. Second, the bus type-based integer programming formulation for heterogeneous fleet problem is more useful in solving real-world problems. Such formulation can reduce the model complexity in terms of the number of decision variables, the number of constraints, and the size of model file. Our benchmark test shows the optimally solvable problem size is increased from 58 trips to 286 trips. At the same time, the fixed costs for each type of buses are considered in the model formulation, which will guide the solver to select as many low priced buses as possible.
Third, we propose a general local search algorithm coupling with simulated annealing to minimize the number of buses and total travel distance. The VRP local search operators, such as one-point move, two-point move, 2-opt move and cross-exchange, are used to improve the route solution gradually. Based on this algorithm, the two-stage metaheuristic solvers for biobjective school bus routing problem or school bus scheduling problem are designed. The test on the RSRB and CSCB instances shows that our algorithm outperforms the existing mixed load heuristic significantly. Compared to the optimal solutions obtained by CPLEX, our algorithm can also provide the optimal numbers of buses for most of the RSRB and CSCB instances. The test on HETERO instances shows that our metaheuristic approach is of obvious advantages for large-size heterogeneous fleet problem in terms of the solution quality and computation efficiency. In addition, the algorithm can be easily integrated with GIS for solving real-world school bus routing and scheduling problem.
For school bus scheduling problem, this paper and most of existing literature focus on the problems with fixed school time windows. However, the benefit of bus scheduling depends largely on the settings of school bell times. For a multi-school system, if the bell times are adjusted properly, the solution for school bus routing and scheduling problems can be improved substantially. Future works might have considered the problems of school bus scheduling and bell time adjustment simultaneously.
Supporting Information S1 File. The benchmark instances files of school bus scheduling problem(SBSP).