A School Bus Routing Heuristic Algorithm Allowing Heterogeneous Fleets and Bus Stop Selection

This paper addresses a school bus routing problem formulated as a capacitated and time-constrained open vehicle routing problem with a heterogeneous fleet and single loads. This problem incorporates several realistic features, such as student eligibility, maximum walking distances, bus stop selection, maximum riding times, different types of buses, multistops, and bus dwell times. A heuristic algorithm based on an iterated local search approach is proposed for this problem. It determines the selection of bus stops from a set of potential stops, the assignment of students to the selected bus stops, and the routes along the selected bus stops. The main objectives are minimizing the number of buses used, the total student walking distance, and the total route journey time. Other aims are balancing route journey times between buses and minimizing the total number of empty seats. A set of 20 real-world problem instances are used to evaluate the performance of the algorithm. Results indicate that the algorithm finds high-quality solutions in very short amounts of computational time.


Introduction
The school bus routing problem (SBRP) is a combinatorial optimization problem that was first investigated over 4 decades ago [25] and has since received a lot of attention within the scientific community (see the section "Related Work"). The problem relates to the organisation of school bus transportation services for students between their home addresses and schools.
School bus operations constitute a challenging task from both a logistics and financial perspective. Typically, an SBRP entails compiling a list of addresses of all students deemed eligible for school transportation, approving potential bus stops reachable by the students, determining the stops to be visited by the buses, assigning students to their respective bus stops, and designing routes that optimize operational efficiency without sacrificing school bus safety and service quality. These objectives are often conflicting in nature, since an improvement in the level of service quality can increase the cost of provision, and vice versa.
In various countries, school bus transportation forms part of the government's administrative mechanism and is funded through taxation. Students who live at least a certain distance from the school they attend are entitled to free or subsidized transport. In Malta, for example, school transport is provided free of charge to all kindergarten, primary, middle, and secondary state school students residing at least 1 km away from their school. Given the large amount of funds being invested in school bus transportation, it is critical for governments to minimize the total cost of these services.
An important objective is to minimize the number of buses, since this reduces their corresponding acquisition and driver employment costs. Another priority is to minimize operational costs by keeping route journey times as short as possible. This also promotes student well-being, particularly for younger children. In Wales, for example, a maximum 45 and 60 min journey time is recommended for primary and secondary school pupils, respectively. Moreover, governments typically issue a policy concerning the maximum distance that students are expected to walk between their homes and designated bus stops. For instance, a walk of 1.6 km is deemed reasonable in Wales.
In this paper, we focus on the single-school SBRP in which routes are constructed for each school separately. Mixed loads, i.e., students from different schools travelling on the same bus at the same time, are not considered, since these are usually not permitted in the locations under study. We dedicate our research to the morning problem, whereby students are picked up from bus stops and dropped off at school. A solution to the afternoon problem, whereby students are picked up from school and dropped off at bus stops, is found by reversing the routes.
The remainder of the paper is organized as follows. The section "Related Work" presents a brief review of related work on the SBRP. The section "Problem Definition" then defines our SBRP, while the section "Algorithm Description" describes our heuristic algorithm for this problem. The section "Computational Experiments" describes the problem instances considered and our computational results. Finally, the section "Concluding Remarks and Suggestions" provides the concluding remarks and some suggestions for future work.

Related Work
The SBRP falls into the class of NP-hard vehicle routing problems (VRPs). According to Laporte et al. [21], VRPs aim to design optimal delivery/collection routes from one or more depots to a finite set of geographically scattered customers under a variety of side constraints. A common objective in VRPs is to minimize the total operating costs of the fleet of vehicles, while typical constraints include maximum capacity restrictions on vehicles [the capacitated VRP (CVRP)] and maximum time/distance restrictions on routes. A taxonomic review of the VRP and its variants is presented by Braekers et al. [8].
The CVRP first appeared over 6 decades ago in the seminal paper of Dantzig and Ramser [14]. In their reviews on SBRPs, both Park and Kim [26] and Ellegood et al. [19] observe that almost all their reviewed publications consider capacity constraints. Here, we take the same assumption but choose to allow for a heterogeneous fleet of buses consisting of several types of buses, with buses of the same type having the same capacity (as in [3,12]). In [19,26], only around 20% and 25%, respectively, of SBRP publications are noted to model a heterogeneous fleet (e.g., [24,29,36,37]). Here, we do the same to extend the homogeneous fleet SBRP that we previously studied in [35] and to shift our research towards more realistic settings. The adaptation of our algorithm in [35] to heterogeneous fleets is the main contribution of this paper.
Another variant of the VRP proposed by Sariklis and Powell [31] is the open VRP (OVRP), in which routes do not start and end at a depot as in the classical VRP, but rather either start or end at a depot. Similar to [2], we model our SBRP as a capacitated and time-constrained OVRP. We do this since governments put out the designed routes to public tender, and thus, buses are kept at bus garages owned by private companies and not at the school. According to [26] and [19], approximately 66% and 35%, respectively, of SBRP publications include maximum riding time constraints (e.g., [11,27]). They also observe that the majority of the publications on school bus routing also deal with the single-school SBRP (e.g., [13,23,33]).
Desrosiers et al. [15] decompose the SBRP into five subproblems. In the first subproblem, data preparation, a network containing the schools, student residences, potential bus stop locations, and bus depots is constructed. Information on the number of students at each residence, the school destination of each student, and the types and number of buses available and their capacities are also specified. The second subproblem, bus stop selection, seeks to select a subset of bus stops from a set of potential bus stops and assign students to these stops. Ellegood et al. [19] remark that less than one-quarter of SBRP research tackles this subproblem (e.g., [6,28,34]). The third subproblem, route generation, is the core of the SBRP. It deals with building initial routes through construction methods such as insertion-based or savings-based heuristics and enhancing them via improvement methods such as metaheuristics and integer programming. The last two subproblems, school bell time adjustment and route scheduling, adjust the schools' opening/closing times to allow buses to service multiple schools and establish chains of routes that can be executed by the same vehicle.
In this paper, we cover the first three subproblems stated above. For the bus stop selection and route generation subproblems, we employ the location-allocation-routing (LAR) strategy (as in [4,5,18]), in which bus stops are first selected, students are assigned to stops, and then, route generation follows. This strategy may lead to sub-optimal routes, since the bus stops selected in the former subproblem may not be best for the latter subproblem [7]. To address this issue, we therefore also include an operator that alters the selection of bus stops and the assignment of students, as discussed in the section "Generation of Alternative Solutions".

Problem Definition
In our problem, we use parameter m w to indicate the maximum distance that students are expected to walk to a bus stop. We also use parameter m e to specify the minimum walking SN Computer Science distance students should live from school to be eligible for school transportation. As in [22,35], our SBRP can be represented via two set of vertices, V 1 and V 2 , and two sets of edges, E 1 and E 2 . The vertex set V 1 consists of one school, v 0 , and n potential bus stops v 1 , v 2 , … , v n , and the edge set E 1 contains all n(n + 1) arcs (u, v) where u, v ∈ V 1 and u ≠ v . Each edge (u, v) in the complete directed graph (V 1 , E 1 ) is weighted by the shortest driving time t(u, v) from vertex u to vertex v. Meanwhile, the vertex set V 2 consists of eligible student addresses, i.e., addresses that are within walking distance greater than m e from the school. Each address w ∈ V 2 has a corresponding number s(w) of students residing at that address and requiring transportation to the school under consideration. The second edge set E 2 is the set It is assumed that the undirected bipartite graph (V 2 , V 1 ⧵{v 0 }, E 2 ) has no isolated vertices. Otherwise, either an address has no bus stop within walking distance m w (and therefore a new bus stop must be added to V 1 ), or a bus stop has no address within walking distance m w and can be removed from V 1 . Moreover, a bus stop v ∈ V 1 ⧵{v 0 } for which there exists an address w ∈ V 2 with the single incident edge {w, v} shall be referred to as a compulsory stop. This is because such a stop v is the only stop within walking distance m w for students living in address w. It must therefore be present in any solution.
A feasible solution to our SBRP is given by a set of routes R = {R 1 , R 2 , …} . An example is shown in Fig. 1. Each route R ∈ R uses one bus which needs to have adequate seating capacity for all the students boarding that route. (The issue of choosing a suitable capacity C(R) for route R will be considered at the end of this section.) This bus successively visits a subset of bus stops and terminates at the school v 0 . The subset of bus stops traversed by at least one route is denoted by This set should cover each address w ∈ V 2 at least once, meaning that students in each address w will have at least one bus stop in V ′ 1 within walking distance m w . Such a covering shall be referred to as a complete covering of V 2 , whereas a covering that does not satisfy this property is an incomplete covering of V 2 . In addition, the total number s(R) of students boarding the bus on route R should not exceed some pre-defined maximum bus capacity C max , and the journey time t(R) of route R should not exceed the maximum journey time m t . These constraints can be expressed as follows: (1) It is important to note that bus stops in V ′ 1 can feature in more than one route in R . For example, there may not be enough spare capacity in a bus to serve all students waiting at a bus stop v ∈ V � 1 . In that case, bus stop v must be visited by more than one bus. Here, a bus stop occurring on more than one route in a solution is called a multistop. In VRP literature, this characteristic is referred to as the allowance of split deliveries [16,17]. It is assumed that each student boarding at a multistop is only permitted to board one specific route serving that stop, since, otherwise, a bus stopping at that multistop may be too full to serve a subsequent stop in its route. Generally, students living in the same address (in particular siblings) are preferably assigned to the same route; however, in practice, one may encounter unavoidable cases where students from the same address are assigned to different routes.
The calculation of the journey time t(R) of route R ∈ R is composed of two components: the total bus travel time and the total bus dwell time. Each dwell time within a route captures the time spent servicing a designated bus stop; i.e., the time spent stopping the bus, opening the doors, boarding the students, and merging back into traffic. In our case, we estimate the dwell time at stop v in route R using the linear function d(v, R) = d 1 + d 2 s(v, R) , where s(v, R) represents the number of boarding students at stop v onto route R, d 2 represents the boarding time per student, and d 1 is a parameter which accounts for the remaining servicing time. Here, d 1 and d 2 are taken to be 15 and 5 s, respectively. Therefore, given a route R = (v 1 , v 2 , … , v l , v 0 ) , the route journey time t(R) is given by where the first two terms give the total bus travel time and the last term gives the total bus dwell time.
The provision of school transportation involves high costs and is sometimes characterized by a shortage of buses. For instance, in September 2018, the Maltese Education Ministry was faced with the issue of having almost 1000 students without transport arrangements due to a lack of space on buses [30]. As previously mentioned, the primary objective of our SBRP is to identify an appropriate subset of bus stops V ′ 1 to minimize the number |R| of routes (buses) included in a solution. In our case, this is achieved by attempting to produce feasible solutions that use the lower bound of ⌈ ∑ w∈V 2 s(w)∕C max ⌉ routes needed to serve all students. Under the assumption of having enough buses of capacity C max to cater for all students, a solution satisfying constraints (1)- (3) and meeting this lower bound of |R| is always guaranteed, since multistops are allowed; however, any one of the routes could potentially violate the maximum riding time constraint (4). Thus, there may be cases where additional routes are needed.
In addition to minimizing the number of routes |R| , in this paper, we also aim to target service efficiency by minimizing the total route journey time in a solution. Whenever multiple feasible solutions have the same minimum total route journey time, we also give preference to the one with the smallest discrepancy between the longest and shortest routes. This encourages equity of the service. Furthermore, we assign students to their closest stop in a selected subset of bus stops. In the future, we would like to incorporate the minimization of student walking distances into the cost function. This is not considered here, but is important for the effectiveness of the service.
At the end of the optimization process, one final task is to assign vehicles to the routes of a solution, such that the total number of empty seats is minimized. This is desirable, because larger buses typically consume higher amounts of fuel. Earlier, we mentioned that an adequate seating capacity is required for each route. In our case, each route is assigned to the smallest pre-defined bus capacity that is greater than or equal to the total number of students boarding that route. We follow this strategy, since, as we said before, governments usually construct the routes before putting them out to tender; hence, each bus type can be assumed to be available in an unlimited number. If buses of each type are limited in number, the problem of assigning buses to routes can be modelled as a minimum-weight maximum-cardinality bipartite matching problem. The bipartition consists of a set of |R| vertices representing the routes in the solution, and a set of B vertices representing all available buses. A route vertex is then linked to a bus vertex if and only if its load is at most the capacity of the given bus. Such an edge is weighted by the difference between the bus capacity and the route load (i.e., the number of empty seats). The Hungarian algorithm, with worst-case complexity O(B 3 ) , can then be used to find a matching that contains as many edges as possible and that has the minimum total weight.

Algorithm Description
Our heuristic algorithm for this problem employs the following overall strategy. To begin, a subset of bus stops representing a complete covering is selected and a nearest-neighbour heuristic is applied to construct an initial solution using a fixed number |R| of routes. As mentioned, |R| is initially taken to be the lower bound stated in the previous section. Note that the initial assignment of stops to routes allows the violation of the maximum riding time constraints (4). A local search routine involving six improvement heuristics is, therefore, invoked on the initial solution to try shortening the routes without changing the current subset of bus stops. After this routine has been completed, a procedure is performed whereby the current subset of selected bus stops is altered and the current solution is repaired. The local search routine is then re-applied. This entire process is repeated until a time limit is reached. If no solution satisfying (1)-(4) is achieved at this limit, then the number |R| of routes is increased by one and the algorithm is restarted.

Construction of Initial Solution
In our approach, the initial subset of bus stops V ′ 1 is selected as follows. First, all compulsory stops are added to V ′ 1 . The non-compulsory stops are then arranged in non-increasing order according to the number of currently uncovered addresses they serve. The stop with the largest such value is then added to V ′ 1 , breaking ties randomly. This ordering and selection procedure is repeated until a complete covering of V 2 is obtained. Each address in V 2 is then assigned to the closest bus stop in V ′

1
. The assignment of addresses to stops determines the number . It may also be the case that some stops have no boarding students, in which case they are removed from V ′ 1 .

Next, each bus stop in V ′
1 is assigned to one of the |R| routes, such that each bus is not overloaded. In our case, this assignment follows a parallel backward implementation of the nearest-neighbour constructive heuristic. To start, |R| empty routes are defined and the remaining capacity c i of each route R i ∈ R is set to C max . The |R| closest stops to the school are then added at the front of the routes, one in each route. Closeness to school is measured by the dwell time at SN Computer Science the stop plus the shortest driving time from the stop to the school. To calculate the dwell time at stop v ∈ V � 1 in route R i , the minimum of c i and s V �

1
(v) is considered as there may be more than c i students boarding at stop v. In this case, a multistop is created, since the remaining s V � 1 (v) − c i students boarding at stop v will need to be assigned to a different route R j . The remaining capacities c i are then updated accordingly. This iterative procedure of determining the closest stop to the most recently added stop in route R i , adding it to the front of the route and updating the remaining capacity c i , is repeated until all stops in V ′ 1 are assigned to a route. On completion, an initial solution R will have been generated and can be evaluated according to the cost function described presently.

Cost Function
Here, a candidate solution R = {R 1 , R 2 , …} is evaluated according to the cost function where This means that if a route R ∈ R satisfies (4), then its journey time is unaltered. On the other hand, if the journey time t(R) of route R exceeds m t , then this journey time is scaled up heavily via a penalty. The addition of the value 1 in the second case of (7) guarantees that two routes both with journey time of at most m t (and, therefore, a cost of at most 2m t ) are always preferred over a single route with a journey time exceeding m t .

Local Search Routine
As mentioned, the intention of our local search routine is to shorten the journey times of routes in a solution R while maintaining the satisfaction of (1)-(3). Note that this local search acts on a solution using a fixed subset of bus stops V ′ 1 . It uses a combination of three intra-route and three interroute local search operators, with the former being applied to any single route R 1 ∈ R and the latter being applied to any pair of routes R 1 , R 2 ∈ R . Without loss of generality, assume that The six operators considered are the following: Cases j = i + 1 and j = i + 2 are special cases of the exchange operator and are thus excluded.
, v j and insert it before stop v k , possibly also inverting the sub-route if this yields a better cost. If k = l 1 + 1 , then the sub-route is inserted before , v j from R 1 and insert it before stop u k in R 2 , possibly also inverting the sub-route if this yields a better cost. If k = l 2 + 1 , then the sub-route is inserted before Here, the value z gives the maximum number of students who can be transferred (hence, decreasing t(R 1 ) as much as possible), such that both occurrences of v i have at least one boarding student and both routes R 1 and R 2 satisfy (3).
The neighbourhood sizes corresponding to the above oper- , respectively. These operators are the same as those used in [22,35]. The exchange, two-opt and crossexchange operators are also used in a similar context in [12], while the generalized Or-opt and Or-exchange are extensions (in cases where i ≠ j ) of operators used in [12,33]. The creating multistops operator was first proposed in [22].
Note that the intra-route operators do not affect the total number of visited stops, the dwell time at each stop, or the total number of boarding students on a route. On the other hand, Or-exchange and cross-exchange moves can lead to a violation of the capacity constraints (3). Such moves are therefore not evaluated. Moreover, the two inter-route operators can result in duplicate stops in the same route, which are removed as follows. Without loss of generality, assume that sub-route v i , … , v j is being transferred from route R 1 to route R 2 and that one stop v h , i ≤ h ≤ j , is already present in R 2 . Then, stop v h is removed from the sub-route and the students boarding this occurrence of v h are all transferred to the occurrence of v h in R 2 .
Our local search routine follows the direction of steepest descent. In each iteration, all moves in the union of the six neighbourhoods are evaluated and the move giving the largest reduction in cost is performed. If multiple moves give the largest reduction in cost, the one which yields the smallest discrepancy between the longest and shortest routes in the solution is performed. Such a breakage of ties aims at balancing the journey times between buses. The local search routine terminates when a solution whose neighbourhoods contain no improving moves is reached.

Generation of Alternative Solutions
As noted, the subset of bus stops V ′ 1 is fixed during our local search routine. However, it may be the case that no solution R for a fixed subset V ′ 1 is feasible with respect to the maximum riding time constraints (4). In this case, altering V ′ 1 and re-running local search may lead to a better solution. For this reason, our algorithm also contains an operator that generates a new subset of bus stops V ′′ 1 , assigns students to these bus stops, and then creates a set of routes that use these stops. We designed four variants of the algorithm, which differ in the way they generate V ′′ 1 . These are the following: I Generate V ′′ 1 from scratch (with no reference to previous iterations); II Generate V ′′ 1 from the subset V ′ 1 used in the previous iteration; III Generate V ′′ 1 from the most recent subset V ′ 1 that yielded a feasible solution with the lowest cost found so far; IV Generate V ′′ 1 via a trade-off between Variants II and III, whereby V ′′ 1 has a 50% chance of being generated according to Variant II and a 50% chance of being generated according to Variant III.
Note that in Variant III, the subset of stops generated in the previous iteration is used if no subset has yet yielded a feasible solution.
In Variant I, the generation of V ′′ 1 follows the same selection strategy as that discussed in the section "Construction of Initial Solution" and new routes are again produced via nearest-neighbour construction. For the remaining variants, the non-compulsory stops in V ′ 1 are identified and a random selection of these is removed. Assuming a total number of non-compulsory stops, in our case, the number of removals is selected according to a Binomial distribution with parameters and 3∕ , so that three stops are removed on average. Upon removal, if we have an incomplete covering of V 2 , then additional stops are added to V ′ 1 . If all addresses not covered by the stops in V ′ 1 are covered by stops that were not originally in V ′

1
, then, at each stage, a stop from the latter set of stops that serves the largest number of uncovered addresses is added, breaking ties randomly. If, on the other hand, some address is uncovered by the stops which were not originally in V ′ 1 , then at least one of the removed stops must be added back. The same selection strategy is applied in this case and the whole procedure is repeated until a new complete covering V ′′ 1 of V 2 is achieved. Each address is then reassigned to the closest stop in V ′′ 1 . As before, stops with no addresses assigned to them are then removed from V ′′ 1 . Having determined a new subset of bus stops, repairs now need to be made to R , so that only bus stops in V ′′ 1 feature in the solution. To do this, all occurrences of stops in (v) students are removed from occurrences of v in R . If this results in an occurrence of v with no boarding students, then this occurrence is removed (v) students can be added, then a new occurrence of v must be added to R . Stops v ∈ V �� 1 ⧵ V � 1 must also be added to the solution. A new stop is inserted in a route having the lowest load, at the position which causes the least increase in the route journey time. If this insertion does not cater for all students boarding that stop, then the procedure is repeated.
Having repaired solution R (or generated completely new routes in the case of Variant I), the local search routine is then re-invoked. This repair-and-improve process is repeated until the time limit is reached.

Computational Experiments
A total of 20 real-world problem instances are considered here, as summarized in Table 1. The problem instances pertaining to the UK and Australia originate from [22] and can be downloaded at [9]. The remainder were generated by us and can be downloaded at [32]. Each problem instance was generated as follows. The location of a school was first identified and a number of random student addresses were selected within the school's catchment area, but more than m e km from the school. The number of students living at each address was generated randomly according to the following distribution: 1, 2, 3, and 4 with probabilities 0.45, 0.4, 0.14, and 0.01, respectively. This distribution approximates the relevant statistics in the locations considered. Potential bus stops were then identified through public records, such that each stop has at least one address within walking distance m w and each address has at least one stop within walking distance m w . The shortest driving times between each bus stop pair and shortest walking distances between each bus SN Computer Science stop and address pair were then determined using the Bing Maps Routes API.
Here, we assume that the maximum journey time m t = 2700 s (45 min) and that vehicles of different capacities are available, depending on the country of the problem instance under study. 1 These capacities were selected by inspecting the fleet of vehicles of several school transport providers in the said locations.

Heuristic Algorithm Results
Our heuristic algorithm was coded in C++ and run on a 3.6 GhZ 8-Core Intel Core i9 processor with 8GB RAM. Variants I to IV were run 25 times on each instance using a time limit of 5 min per run. Overall, we found that feasible solutions using the lower bound of ⌈ ∑ w∈V 2 s(w)∕C max ⌉ routes were achieved in nineteen of the 20 instances in all runs. The only instance which required one additional route was the Bridgend instance. This is probably because the pairwise travel times have 25th, 50th, and 75th percentiles approximately equal to 9.4, 14.4, and 18.9 min, respectively. Since these percentiles are relatively high, there is a large chance of violating the maximum riding time constraints (4). Table 2 displays summary statistics on the results achieved by our algorithm for heterogeneous fleets of buses. Columns 3 to 6 display the average number of iterations performed by each algorithm variant. For the Suffolk instance, approximately 99.9%, 99.9%, 69%, and 81.3% of the iterations performed in Variants I to IV, respectively, produced infeasible solutions, on average. For Cardiff, these percentages were 12.6%, 45.6%, 25.5%, and 33.3%, while those for Bridgend were 13.8%, 6%, 0.1%, and 3.7%. Some infeasible solutions were also produced in Variants II and IV for Porthcawl (0.04% and 0.02%, respectively) and Variant I for Brisbane (0.005%, 1 iteration only). Constraints (1)-(4) were satisfied in all runs for all remaining 15 instances.
According to Table 2, Variant II performs the highest average number of iterations for all instances except and Edinburgh-2. As expected, Variant I performs the lowest average number of iterations for all instances except Porthcawl. This is because this variant does not use information from previous iterations when altering the current subset of bus stops. Hence, local search takes longer in each iteration as it operates on a newly constructed set of  Fig. 2. In each boxplot, the horizontal axis represents the total journey time, the endpoints of the whiskers represent the minimum and maximum of all journey times, the red box displays the variability between the 25 and 50th percentiles, and the grey box displays the variability between the 50 and 75th percentiles. Figure 2 indicates that the performance varies significantly across variants for all instances except Mġarr. For the latter instance, all variants provided the same total journey time in all runs. Variants III and IV are seen to perform best for the majority of the instances.
To demonstrate that the results of the four algorithm variants are statistically significantly different, Kruskal-Wallis and post-hoc Bonferroni-adjusted pairwise comparison tests were performed. As expected, significant differences between the variants (all p values < 0.001 ) were revealed for all instances except . The pairwise comparison tests indicated that, for 16 instances, the total journey times of Variants I and II are significantly different at the 0.05 level of significance than those of Variants III and IV. Ten of these instances also saw a significant difference between Variants I and II. Two instances (Porthcawl and Edinburgh-2) saw significant differences between Variant III and all other variants, whereas the Cardiff instance saw significant differences between Variant II and all other variants.
Moving on to Table 3, Columns 5 to 8 display each instance's best total journey time for Variants I to IV, respectively. Each best total journey time corresponds to the endpoint of the left whisker in the boxplot (refer to Fig. 2). Each instance's best reported result across all variants is displayed in bold and the number of runs giving that result is shown in brackets. According to the table, Variants I to IV produce the best reported results in 3, 4, 15, and 12 instances, respectively. Variants I and II seem to be most appropriate for small-sized instances, Variant IV for small-to-moderatelysized instances, while Variant III is more effective for large instances. Furthermore, Variants I to IV produced best total journey times that are at most 13.93%, 10.54%, 4.28%, and 1.31%, respectively, worse than the best reported results. Our best reported result for eleven instances was achieved in only one run, whereas multiple runs reached the best reported result for the other nine instances. Some or all multiple runs for all the latter instances except Porthcawl and Edinburgh-2 also have different corresponding subsets of bus stops. The total number of alternative subsets of bus stops is given in Column 4.  Table 3 gives the bus capacities assigned to the routes of the best reported result for each instance. The superscript displayed after each capacity indicates the required number of buses of that capacity. It must be pointed out that the capacity assignment of different solutions yielding the best reported result is not necessarily the same. For example, the two solutions for the Senglea instance have different assignments; one uses one 14-passenger bus and five 53-passenger buses, while the other uses one 18-passenger bus, one 44-passenger bus, and four 53-passenger buses. In such a case, the assignment which provides the lowest number of empty seats is presented. This number is displayed in Column 3 of Table 3.
Recall that our algorithm follows the steepest descent (SD) direction, meaning that the best improving move from all neighbourhoods is executed in each iteration. To evaluate whether a change in this strategy would significantly affect the quality of the solutions obtained, the algorithm was altered, such that it follows a variable neighbourhood descent (VND) direction. Specifically, in each iteration, the first improving move from a randomly selected neighbourhood out of the six considered is executed. The process continues until all neighbourhoods do not contain improving moves. As expected, this alteration led to an increase in the average number of iterations performed by each algorithm variant for the majority of the instances. In fact, 18 of the 20 instances saw a percentage increase ranging between 0.44 and 124.88% for all four variants, implying that VND iterations take less time on average. The Porthcawl and Edinburgh-2 instances again stand out from the rest, with the former having percentage increases in the range [− 0.30%, 0.04%] and the latter in the range [− 1.45%, − 1.24%]. When the best total route journey times of the SD and VND strategies were compared, it was observed that these were the  Table 3 Best total route journey times (in min) achieved by our heuristic algorithm across 25 runs Column 9 (Inc.) gives the total route journey time of the best incumbent solution achieved by Gurobi. Column 10 displays the percentage improvement between the best reported result from our heuristic algorithm and the result in Column 9. The last column gives the relative gap between the result in Column 9 and the best-known lower bound on the optimal total route journey time Location  (19) 67.12 (5) 67.12 (25) 67.12 (25)  same for eleven instances. For five instances (Pembroke, Canberra, Valletta, Birkirkara, and Brisbane), the VND's best solution was worse by 0.02-1.86%. On the contrary, the Victoria, , Edinburgh-1, and Adelaide instances saw superior best solutions from the VND strategy by 0.04-1.21%. In addition, a sign test was performed for each instance to determine whether the results of the two strategies are statistically significantly different at the 0.05 level. Significant differences were found for seven instances, with registering better results from VND, and Qrendi, Pembroke, Canberra, , Birkirkara, and Brisbane registering better results from SD. This suggests that, in our case, SD overall performs better than VND.

Mixed Integer Programming Results
An attempt was also made to improve each instance's best reported result from our heuristic algorithm. For this purpose, a mixed integer programming (MIP) model was formulated for our SBRP, as shown in the appendix. This model was executed using Gurobi 9.1 with a run time limit of 2 h per instance.
For each instance, the best reported solution from the heuristic algorithm was provided as an initial solution. For instances with multiple best solutions, one solution was chosen randomly. Warm-starting the MIP in this way helps to speed up the convergence of branch-and-cut, which is the method used by Gurobi.
The MIP results are presented in Columns 9 to 11 of Table 3. Column 9 gives the total route journey time of the best incumbent solution found. Column 10 displays the percentage improvement between the best reported result from our heuristic algorithm and the best incumbent result. Finally, Column 11 gives the relative MIP optimality gap between the best incumbent result and the best-known lower bound on the optimal total route journey time. We see that Gurobi was able to make improvements for seven instances and the percentage improvements of these range between 0.02 and 5.55%. On the other hand, Gurobi was not able to improve the heuristic algorithm solution of nine instances within the time limit. For the remaining four instances, the solver could not provide any results before the time limit was reached.

Comparison to Existing Benchmarks
Further experimentation on the performance of our heuristic algorithm was carried out by applying it to a set of 100 instances generated by Sales et al. [29]. In these instances, the sizes of V 1 ⧵ {v 0 } and V 2 range from 25 to 250 bus stops and from 500 to 5250 addresses (or students, as each address is assumed to correspond to one student). Six assumptions in [29] that differ from the ones we employ in this paper are the following: (i) routes are defined as cycles that start and end at the school, (ii) the heterogeneous fleet is limited in size, (iii) fixed costs are associated with bus usage which vary according to the bus capacities, (iv) driving distance is taken into account instead of driving time (hence, no bus dwell time is incorporated), (v) no maximum limit on the routes' length is considered, and (vi) students are assigned to the stop, within maximum walking distance, that is closest to the school. Moreover, the objective in [29] is to minimize the sum of the total distance travelled by all buses and the fixed costs associated with the buses used.
In [29], the instances were tackled via a memetic algorithm that we describe briefly. Initially, each student is assigned to a stop and the stops with at least one student are identified. An initial population of 100 solutions is constructed by randomly selecting the initial stop for each solution and then, at each step, appending a random stop from the three closest stops to the previously visited stop. Buses are then selected for each solution using the following strategy. If there are k buses with capacities C 1 , C 2 , … , C k and fixed costs F 1 , F 2 , … , F k , then the bus i ∈ {1, 2, … , k} which minimizes (C i − D i )F i is selected, where D i is the highest accumulated number of students (following the order of the stops in the solution) not exceeding C i . Subsequently, two parent solutions are selected through the binary tournament method, and the edge recombination crossover operator is applied to these solutions to create an offspring. Up to 50 attempts are made to improve the offspring via the lambda interchange local search procedure with = 2 . If the offspring is better than the worse parent, then it replaces the latter in the population; if not, the offspring still has a 2% chance of replacing the worse parent. This process is repeated for 5000 iterations.
The following are the changes that were made to our heuristic algorithm to test it on the above-mentioned instances: • The cost function (6) was changed as in [29]. In particular, we incorporated fixed bus costs and considered driving distances instead of driving times. The bus dwell times and maximum riding times were removed, and routes were also defined as cycles starting and ending at the school, rather than paths starting at a bus stop and ending at the school. • Students were assigned to the stop, in a selected subset of bus stops, that is closest to the school, rather than to the one that is closest to their residence. • Given that the fleet size is limited, to verify that an inter-route operator move is feasible, the capacities of the buses performing the two routes as well as the capacities of the unused buses were checked to see whether there are two sufficient capacities for the new loads; • The creating multistops operator was considered for any pair of distinct routes R 1 and R 2 , attempting each possible transfer amount z between 1 and the difference between the maximum of the capacity of the bus performing R 2 and the highest unused bus capacity, and the load of R 2 ; • In Variants II to IV, whenever the total number of non-compulsory stops was positive and at most three, the parameters of the Binomially distributed number of removals were changed to and 1∕ , and the solution alteration process was skipped whenever = 0. The computational time allowed for each run of our algorithm was selected as follows. First, instances for which an average computational time of fewer than 10 s was reported in [29] were run for 10 s. For the remaining instances, the average computational time reported in [29] was rounded up to the nearest 30 s and set as a time limit. This was done to have similar running times for comparison purposes. We observe in Table 4 that the results obtained by our algorithm are very close to those obtained by the memetic algorithm, and sometimes better. The relative gaps between the results range between − 6.25 and 4.69%. We have found better results for 68 of the 100 instances, with the largest improvements being seen for Instances 61, 12, and 22. Our algorithm performs the worst in comparison with the memetic algorithm for Instances 40, 30, and 20. Table 4 shows an overall pattern in the relative gaps achieved for different maximum walking distances m w . One can note that, in most cases, the relative gaps increase with an increase in m w . In fact, our algorithm always performs better for m w ∈ {5, 10} . For m w ∈ {15, 20, 25} , our algorithm performs better for 17, 7, and 4 instances, respectively. This trend is likely to have occurred due to the different order of the bus stop selection and student assignment strategies. While our algorithm first selects bus stops and then assigns students to the stop closest to the school, the memetic algorithm performs the student assignments before the selection of the bus stops. Hence, for large m w , as Sales et al. have observed, the memetic algorithm "can concentrate students at stops closer to the school", thus reducing the total distance travelled by all buses.

Concluding Remarks and Suggestions
This paper has addressed a real-world SBRP that incorporates heterogeneous fleets and bus stop selection. This problem includes several important features, such as student eligibility, maximum walking distances, maximum riding times, different types of buses, multistops (multiple buses visiting a single bus stop), and bus dwell times. A heuristic algorithm has been developed which encompasses the first three subproblems of the SBRP, as defined in Desrosiers et al. [15].
Our heuristic algorithm has been tested on a variety of real-world problem instances from Malta, the UK, and Australia, with sizes upwards of 1800 potential bus stops and 750 students. For all instances, the algorithm has successfully found high-quality solutions in short amounts of computational time. A merit of our algorithm is that it can sometimes provide multiple subsets of bus stops that yield the same best total journey time. Indeed, here, we have determined alternative subsets for seven instances. From these alternative subsets, the most appropriate one can then be identified based on factors, such as bus depot locations, bus stop accessibility, and average student walking distance.
We have attempted to improve the best reported results of our heuristic algorithm through an MIP. An improvement, ranging between 1 and 197 s, has been found for seven of the 20 best results. Higher improvements could potentially be obtained by allowing a longer time limit for each instance.
We have also applied our heuristic algorithm to 100 instances generated by Sales et al. [29], with sizes up to 250 potential bus stops and 5250 students. We have observed that our algorithm compares well with the memetic algorithm presented in [29] and also provides better results for 68 instances.
In the future, we would like to handle uncertainties in the bus travel times caused by factors such as weather conditions and traffic congestion (e.g., [1,10,38]). Due to these uncertainties, solutions may no longer remain the best or even feasible, since, for example, students may arrive late for school. It is of interest to see how the solutions will change when we apply a robust optimization and/or a stochastic programming formulation to our SBRP. Another interesting future development is the incorporation of time windows for arrivals at schools and multi-tripping. In the latter, several routes, possibly pertaining to different schools, are linked, so that buses can perform multiple routes successively. Relative gaps between the two results are also presented for each instance The following MIP model produces solutions consisting of cycles that start and end at the school. The arc from the school to the first bus stop in each route is then excluded. This is made possible by assuming that the driving time from the school to any stop is zero. The decision variables of our model are as follows. Binary variable x uvR indicates whether route R ∈ R travels from u to v, where u, v ∈ V 1 and u ≠ v . Binary variable y vR indicates whether route R ∈ R visits v ∈ V 1 . In addition, binary variable z wv indicates whether students in address w ∈ V 2 walk to stop v ∈ V 1 ⧵{v 0 } . Integer variable s vR ∈ {0, 1, … , C max } gives the number of students boarding route R ∈ R at stop v ∈ V 1 ⧵ {v 0 } . Moreover, integer variable l vR ∈ {0, 1, … , C max } gives the total load of route R ∈ R just after visiting stop v ∈ V 1 ⧵{v 0 } . Finally, the continuous variable t R ∈ [0, m t ] specifies the total journey time of route R ∈ R . The MIP formulation is as follows: y vR ≤ s vR ∀v ∈ V 1 ⧵ {v 0 }, R ∈ R, Here, the objective function (8) minimizes the total journey time of all routes. Constraints (9)-(11) relate to stop and school visits. Specifically, (9) and (10) ensure that if route R ∈ R visits v ∈ V 1 , then route R should enter and leave v exactly once. Next, (11) forces each route R ∈ R to visit school v 0 whenever it visits at least one stop v ∈ V 1 ⧵{v 0 } . Constraints (12)- (14) relate to student walks and pickups: (12) ensures that students living in each address w ∈ V 2 walk to exactly one stop within walking distance m w ; (13) ensures that no student walks to an unvisited stop; while (14) guarantees that the total number of students boarding at stop v ∈ V 1 ⧵ {v 0 } is equal to the total number of students walking to that stop. Constraints (15)-(16) relate to student boardings. These force the number of students boarding route R ∈ R at stop v ∈ V 1 ⧵ {v 0 } to be zero if route R does not visit stop v. If route R visits stop v, then (15) also updates the lower bound on the number of boarding students to one. In addition, Constraints (17)-(18) relate to route loads and also serve as subtour elimination constraints as proposed in [20]. Note that l v 0 R = 0 ∀R ∈ R . These guarantee that no route contains a subtour disconnected from school v 0 and that each route load increases in accordance with the number of students boarding the bus on that route. In fact, if route R ∈ R goes from u ∈ V 1 to stop v ∈ V 1 ⧵ {u, v 0 } , then the load of route R immediately after visiting stop v is set equal to the load of route R just after visiting u plus the number of students boarding route R at stop v. Finally, (19) calculates the total journey time of each route R ∈ R.