The cumulative school bus routing problem: Polynomial‐size formulations

This article introduces the cumulative school bus routing problem, which concerns the transport of students from a school using a fleet of identical buses. The objective of the problem is to select a drop‐off point for each student among potential locations within a certain walking distance and to generate routes such that the sum of arrival times of all students from their school to their homes is minimized. The article describes six polynomial‐size mixed integer linear programming formulations based on original and auxiliary graphs, and the formulations are numerically compared on real instances. The article reports the results of computational experiments performed to evaluate the performance of the proposed models.


INTRODUCTION
School buses are the most preferred and frequently used transport system of commuting for students.There are several advantages to transporting students by bus, including safety, maintaining punctuality, promoting student activity, as well as environmental benefits associated with reducing traffic and pollution.Several operational decisions are associated with routing school buses, such as travel routes, stop selections, and stop numbers, collectively referred to as school bus routing problems.Following what appears to be the first study by [23] on school bus routing, numerous studies have appeared to address various issues.
Of particular interest is the study by [5] which lists efficiency, effectiveness, and equity as important criteria in school bus routing problems.Efficiency takes into account service costs, while effectiveness measures a level of service for all eligible students.Equity relates to fairness of the service among the students, which is a crucial but often neglected issue in school bus routing problems [27].
In this article, we describe a variant of the school bus routing problem, namely the cumulative school bus routing problem (CSBRP), with a focus on service times for students.The problem consists of using a fleet of identical school buses to plan for the daily transport of students from one school to their homes.The problem involves selecting a drop-off point for each student that is within a certain walking distance from the student's home, as well as generating routes that start from the school and visit the selected drop-off points to deliver the students to their homes.Selecting the stop points and generating the routes should be such that the sum of arrival times for the students at their homes is minimized.We consider a fixed time for the bus departure from school.The objective of the problem is to minimize the total (i.e., the sum of) arrival times of the students at their homes when travelling from their school.This is equivalent to minimizing the average arrival time.For a fixed departure time from school, the arrival time for each student is determined by the time spent on a school bus plus the walking time from a drop-off point to the student's home.The objective of the CSBRP differs from the classical sum of route lengths, which is intended to FIGURE 1 CSBRP characteristics and details of a route visiting bus stop b 1 and student homes s 3 and s 4 .Students s 1 and s 2 are dropped off at b 1 and walk home from there.Similarly, student s 5 is dropped of at s 4 and walks home from there.minimize the operator's costs.According to [7], minimizing the total route distance may not properly reflect the need for equity and fairness.
In order to illustrate the characteristics of the CSBRP, we find it useful to provide a small numerical example.Figure 1 shows an example of CSBRP characteristics and elements related to a single route in a solution.
In Figure 1, the solid arcs show a bus route which visits bus stop b 1 followed by student homes s 3 and s 4 .Each dashed arc represent walking from a drop-off point to home.The number at each solid arc shows the travel time for the bus along that arc, and the number at each dashed arc shows the walking time along that arc.For each student, the number in parenthesis is the arrival time at the home node, assuming that the route begins at time 0 at the school.
For students s 3 and s 4 , each arrival time is given by the travel time for the bus until the home node.For students s 1 , s 2 , and s 5 , each arrival time is given by the bus travel time until the drop-off point plus the walking time from the drop-off point to the home node.
The example illustrates the three possibilities for each student with respect to drop-off point.Specifically, a student can be dropped off at a bus stop (as for s 1 and s 2 ), at the home node (as for s 3 and s 4 ), or at another student's home (as for s 5 ).Both students s 1 and s 2 get off the bus at b 1 , student s 3 gets off the bus and the home node, and both students s 4 and s 5 get off the bus at s 4 .
The sum of the five arrival times is 49, which is the value of our proposed objective function.
We would like to emphasize that although this example shows only a single route, the CSBRP generally involves determining multiple routes for different vehicles, all beginning at time 0 at the school.The calculation of arrival times is done for each route as shown, and the resulting objective function value is obtained by summing the arrival times over all students.
The goal of this article is to provide mathematical formulations for the CSBRP and solve the problem with a real-life data set.To this end, we propose and computationally compare six polynomial-size mixed-integer linear programming formulations for the CSBRP, based on different definitions of graphs and decision variables.
The remainder of this article is organized as follows.A review of the relevant literature presented in Section 2. In Section 3, the six mathematical formulations are provided.Computational analysis, which consists of data preparation, results, and sensitivity analysis, is presented in Section 4. In Section 5, we address and motivate our choice of unit-demand formulations and describe the model changes if we had used general demands.Conclusions with suggestions for future research are given in Section 6.

LITERATURE REVIEW
A comprehensive literature review of the SBRP until 2009 is conducted by [27], and the review with contemporary research is updated by [10], showing that this problem still attracts the attention of many researchers from the first time that it was studied by [23].
The SBRP is a special case of pickup and delivery problems, as the students are picked up from different locations and delivered to one location.It should be noted that pickup and delivery of people, referred to as dial-a-ride problems (DARP), is different from pickup and delivery of goods since passenger convenience must be considered as well as operational costs.Literature reviews of the various models and methodologies associated with the DARP can be found in [8,15].
The SBRP is also similar to the m-Ring Star Problem (m-RSP), introduced by [2], which has applications in telecommunication problems.The goal of the m-RSP is to provide m tours by either visiting customers or transition nodes, so that the total route cost, together with assignment cost for customers who are not visited on any tours, can be minimized.However, in the m-RSP, the customers without an assigned node in the tour can be assigned to any nearest node in the tours, while in the SBRP, students' potential drop-off points are restricted to the stops that are located within students' allowed walking distance.Moreover, our proposed objective is different from the m-RSP.
According to the study by [27], the main sub-problems of SBRPs are described as follows: bus stop selection, bus route generation, bus route schedules, school bell time adjustment, and strategic transport policy.This article focuses on a combination of two sub-problems: selecting stops and planning efficient routes for a fleet of vehicles.The combination of these two sub-problems is also categorized as location-routing problems in the literature (see, e.g., [22]).
In summary, according to [27], SBRPs are studied with different objectives and constraints, and sometimes they are considered together or interchangeably.The common objectives are to minimize the number of buses, the total travelled distance by all vehicles, or the total walking distance of all students.Some of the studied constraints are limitations on maximum allowed walking distance of students, maximum pick-ups at each stop, or maximum operating time for each vehicle.As an example, [3] proposed an integer linear programming (ILP) model to solve the real-world single school bus routing problem for transporting students of an elementary school in central Ankara, Turkey.They considered a capacity constraint for vehicles and a maximum travel distance constraint for vehicles, with the objective of minimizing the bus operating cost.Other studies with the objective of minimizing the total route costs are [14,18,33].
There are some objectives and constraints in SBRPs that are in conflict.For instance, reducing the total distances travelled by buses results in the students walking longer distances.As such, some of the articles address these conflicting terms in their objective like in the studies by [5,6,29,31].
Five mathematical models were provided by [30], which were the three-index variables formulation, the set partitioning formulation based on the Dantzig-Wolfe decomposition, and the next three models enriching the set partitioning by adding additional inequalities and defining compact versions of the set partitioning formulation.The objective of their proposed model was to minimize the total route distance by vehicles and the total walking distance of students.For effective load balance between routes, a given maximum number of stops per vehicle and a given maximum distance and time that students spent in each vehicle were considered as constraints.
In the study by [26], the objective is to minimize the maximum bus ride length in order to increase equity between passengers.However, in relation to school bus routing, it should be noted that minimizing the maximum ride length only minimizes the maximum arrival time at a student, whereas the CSBRP, which is the objective of our proposed model, also considers the arrival times at all intermediate students, which can vary considerably if many students are served early or late on their routes [24].
However, there is little research on minimizing the total arrival time for all students in SBRPs.The arrival time is considered separately for each student, and it differs from the total travel times of the vehicles.The accumulated arrival times include the time students spend on school buses from school to drop-off points, as well as the walking time of the students from the drop-off points to their homes.In other words, the accumulated arrival times is the sum of the arrival times of all students from a school to their homes.This type of consideration is also referred to as fairness among all students to get to their home after school since it involves the arrival times of all students.
One of the first articles that bring the idea of cumulative travel distance was introduced by [11] in the delivery man problem, where an objective was to minimize the total arrival times to the customers.This type of problem arises when priority is given to customer satisfaction.A new formulation for the cumulative capacitated vehicle routing problem (CCVRP) was defined by [24], which is a generalization of the travelling repairman problem by including the capacity for a fleet of homogeneous vehicles, and in the literature, it is also known as a minimum latency problem.In another study, [20] proposed the set partitioning formulation for the CCVRP, and they solved the problem by a branch-and-cut-and-price algorithm.
Minimizing the cumulative arrival times is also critical in delivering humanitarian aid when disasters strike, in order to save lives and ensure that essential supplies arrive fast [7].
In this article, we describe and numerically compare six compact mathematical models for the CSBRP, which are different in graph definition and decision variables, in the same spirit as some similar studies in the literature that compare compact models.For instance, the classification and comparison of polynomial formulations for the asymmetric travelling salesman problem have been studied by [25,32].
In the study by [34], the authors provided four compact formulations for the generalized travelling salesman problem (GTSP) with time windows, which has an application in last-mile delivery problems.In their study, a set of nodes are partitioned into a number of clusters, which represent the possible delivery locations associated with each customer, and the objective is to find the minimum travel cost by a single vehicle, such that each cluster is visited exactly once.They compared four models based on results obtained by linear relaxation and the branch-and-bound scheme which was implemented in CPLEX.Unlike their research, our proposed model examines another objective of SBRP.Moreover, in our proposed model, the drop-off points for the students are not totally separated, since some of the stops can be shared among several students within a defined walking distance.
Finally, we compare our proposed CSBRP to the capacitated vehicle routing problem (CVRP), the capacitated open vehicle routing problem (COVRP), and the CCVRP.The COVRP is a variation of the CVRP in which a fleet of vehicles is assigned to deliver goods or services to a set of passengers or locations.In the COVRP, unlike the traditional CVRP, each route starts at the depot and ends at any passenger, rather than requiring each route to return to the depot after servicing the last passenger on the route.In this case, the cost of traversing an edge between a last visited passenger and the depot is not included in the objective.The CSBRP, which is studied in this article, may be considered as a variant of the CVRP, as the beginning and ending points of the routes are at schools.However, as the objective of the CSBRP is to minimize the sum of arrival times for all students, the time that a school bus will need to return to a school after the last visited drop-off point is not included in the objective, which is similar to the COVRP.Moreover, the cumulative nature of the objective function makes the CSBRP different from both the CVRP and the COVRP, whereas the CCVRP, on the other hand, also has a cumulative measure in the objective function, as mentioned later in this section.However, the multi-modality as given by both bus transport time and walking time in the objective function of the CSBRP, as illustrated in Figure 1, is not considered in any of the problems CVRP, COVRP, or the CCVRP.As such, the introduced CSBRP share some characteristics with the CVRP, the COVRP, and the CCVRP, while still being different from all of them.

MATHEMATICAL MODELING
This section defines and presents six polynomial-size mathematical models for the CSBRP, with two based on the original and four on an auxiliary graph.Those defined on the original graph use three-index variables, using, for subtour elimination, either Miller-Tucker-Zemlin (MTZ) type constraints [21], or those based on single-commodity flow in line with Gavish and Graves (GG) [13].The formulations based on the auxiliary graph use either two-or three-index formulations, again using either MTZ or GG-type subtour elimination constraints.The models on the original graph are explained in Section 3.1, and the models on the auxiliary graph are explained in Section 3.2.

The CSBRP models on the original graph
The CSBRP is defined on a directed graph G = (N, A), where N is the set of nodes, which contains a school node {0}, a set S of student nodes, and a set B of bus stops, and A = {(i, j) ∶ i, j ∈ N, i ≠ j} is the set of arcs.A fleet of identical school buses, the index set of which is denoted by K, each of capacity of Q, start from the school and drop off all students picked up from the school.The travel time between any two nodes is denoted by T ij for any (i, j) ∈ A. For each student s ∈ S, the subset U s ⊂ N shows a set of potential stops reachable by student s, and are the potential drop-off points (the student's home location, other students' home locations, and bus stops) within an allowable walking distance W Max from the student's home.We denote W si as the walking time of student s ∈ S between node i ∈ U s and the home address of student s.
When W Max = 0, students must be dropped off at their home addresses, and when W Max is greater than zero, students may be dropped off at locations other than their home addresses.In this case, some potential stop points will likely be shared among a number of students.For instance, when stop k is located within the walking distance of both students s and t, both subset U s and U t will contain stop k.If a drop-off stop is shared between a number of students, each should have the same opportunity to be dropped off there at any time by school buses.By this assumption, it is possible that one stop is visited by a number of vehicles at different times, just as it is possible to deliver multiple students by the same vehicle in one visit to a stop point.This possibility is also referred to as multiple visits of one node.
In the following subsections, we are going to provide self-contained models for the CSBRP.For each model, we first define the decision variables and then the related formulation.

Model 1: Three-index formulation on the original graph based on MTZ
The following is the definition of decision variables that are used for the three-index formulation on the original graph, and the subtour eliminations are based on MTZ constraints.

Decision variables
x k ij is a binary variable, which is equal to 1, if arc (i, j) ∈ A is traversed by vehicle k; otherwise, it is zero.
z k sj is a binary variable which is equal to 1 when student s ∈ S is dropped off at stop j ∈ U s by vehicle k; otherwise, it is zero.
is a non-negative value that shows the arrival time at node i of vehicle k, if vehicle k visits node i; otherwise, it is zero.
The objective of the model is: which calculates the arrival time of each student separately.As explained earlier, it is also possible that several students are dropped off at one node by the same vehicle or by different vehicles at different times.Therefore, we cannot uniquely determine and use the time of visit for each node.
In order to make the objective function (1) linear, we define a new decision variable t s which is a nonnegative value and shows the arrival time at the drop-off point for student s.
The objective function ( 2) is to minimize the sum of arrival times for all students that contains the time that all students spend on the vehicle up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.Constraint (3) shows that each student is dropped off at only one of the stops in the defined subset.Constraint (4) shows the arrival time of student s to node i by vehicle k, where M is a sufficiently large value.Constraint (5) shows the number of vehicles that depart from the school.Constraint (6) ensures that if vehicle k visits a stop, it leaves the visited stop.Constraint (7) shows the flow of each vehicle at each node i. Constraint (8) avoids the assignment of students to a non-visited stop.Constraints (9) ensures that the maximum number of onboard students on each vehicle cannot exceed the vehicle capacity.Constraint (10) shows the arrival time consistency and imposes the time connectivity for each vehicle, and also eliminates subtours in the same way as shown in [9].Constraints ( 11)- (15) show the domains of the decision variables.
Constraint (10) can be strengthened as below, as shown in [9] for the distance-constrained VRP: where T is a sufficiently large value, for which a bound is given below:

is the ith largest travel time between any pair of nodes in N.
Proof.We seek a small value for T such that it does not cut off any optimal solutions.If both x k ij and x k ji are zero, then Constraint (16) reads As a result, we are looking for the largest value for the arrival time, or more precisely a relatively small value that is larger than any arrival times within the routes.The largest arrival time within routes is also the same as the arrival time of the last student who is dropped off by the school bus.Considering the school buses' capacity, the maximum number of students on the school buses is limited to Q.If there are Q students on board, the arrival time of the last student on a route is less than or equal to the sum of the Q largest travel times between any pair of arcs.▪ 3.1.2Model 2: Three-index formulation on the original graph based on GG In this formulation, we are going to show the three-index formulation on the original graph based on GG inequalities.

Decision variables
x k ij is a binary variable, which is equal to 1, if arc (i, j) ∈ A is traversed by vehicle k; otherwise, it is zero.
z k sj is a binary variable, which is equal to 1 when student s ∈ S is dropped off at stop j ∈ U s by vehicle k; otherwise, it is zero.
is a non-negative value that represents the number of students who traverse arc (i, j) ∈ A by vehicle k.If arc (i, j) is not traversed by vehicle k, then the value is zero. s.t.: The objective function (17) minimizes the sum of arrival times for all students; this contains the time that all students spend on the vehicle up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.Constraint (18) shows that the students are dropped off only at one of the stops in their subset.Constraint (19) shows the number of vehicles that depart from the school.Constraint (20) ensures that if vehicle k visits a stop, then it leaves the stop.Constraint (21) shows the flow of each vehicle at each node i. Constraint (22) avoids to assign students to a non-visited stop.Constraint (23) ensures that vehicle k must drop off a student in node i if the vehicle visits node i. Constraint (24) ensures that at each visited stop, at least one student is dropped off.Constraints ( 25)- (27) show the maximum number of on board students on each vehicle.Constraints ( 28)- (31) show the domains of the decision variables.
Note that the objectives of function (2) and function (17) are the same, but they calculate the sum of the arrival times of the students in different ways.The objective function ( 2) is keeping track of the arrival time of the students separately, while the objective function ( 17) is calculating the total arrival time of the students by multiplying the number of onboard students who are passing the same arc by the travel time of the arc.

The CSBRP models on an auxiliary graph
As explained earlier, we assume that multi-visits of the nodes are allowed in our models.One of the ways to model this condition is to decompose the graph into clusters with a separate cluster for each student, and impose the requirement that each node can be visited at most once in each cluster.To decompose the graph, we first duplicate nodes to obtain a graph where any node, except for the school node, can serve as a drop-off node for only one specific student.We call this new graph the auxiliary graph, where the same physical location may be represented by several nodes in different clusters.This is the case whenever a node in N has been duplicated.We then partition the nodes in the auxiliary graph into clusters, similar to the generalized vehicle routing problem (GVRP) as introduced by [12], which is a generalization of the capacitated VRP.The GVRP is also a generalization of the GTSP, for which ILP formulations appear in [12,19].The GVRP takes as input a graph partitioned into clusters, and determines the routes that start and end at the school, visit exactly one node in each cluster, and yield the minimum total travel cost see, for example, [1,4].
The auxiliary graph G H = (N H , A H ) is constructed as follows.For each student j ∈ S, we create the cluster C j as a copy of the node set U j ⊆ N of potential stops for student j and define N H = ∪ j∈S∪{0} C j .Therefore, the node set N H is partitioned into (nonempty and disjoint) |S| + 1 clusters, in which cluster C j consists of the home address of student j ∈ S and copies of other nodes from N that are within the allowed walking distance (W Max ) from student j's home address.As such, a distinct node in N H is used for each combination of student and drop-off location.For each student j ∈ S and each node i ∈ C j , we let W i denote the walking time between the home address of student j and node i.In addition, we define W 0 = 0. Note that W i is the walking time, and it differs from W Max , which is the maximum walking distance that the students are allowed to walk, and W Max is used in finding the potential stops for the students in the models.Since two nodes i and k in different clusters may represent the same physical location but for two different students, W i and W k for such a pair of nodes may very well take different values, representing the two students' individual walking times from the shared physical drop-off location to their respective homes.Generally, the same location, that is, node in N, will be represented by a separate node in each cluster which represents a student that can be dropped off at that location.
The arc set u ≠ v} contains arcs that only connect nodes between different clusters, with T ij denoting the travel time on arc (i, j) ∈ A H and T ii ′ = 0 for node i and its copy i ′ in a different cluster as they share the same physical location.
Figure 2 shows a schematic view of the potential stops for the students in both the original graph and the auxiliary graph, where there are seven nodes, that is, the school (0) and three students (1, 2, and 3) displayed by circles and three bus stops (4, 5 and 6) displayed by squares.Figure 2A shows that there are four subsets U j for j = {0, 1, 2, 3} in the original graph, one for the school and the rest are the potential stops for the students.The dashed circle around each student shows the maximum allowed walking distance, and the nodes within each subset indicate the potential drop-off points for the corresponding student.For instance, the drop-off points for student 1 are the home address of student 1 and bus stop 4. Both students 2 and 3 can be dropped off at bus stop 5, and at the home address of student 2 or 3.In addition, student 2 can be dropped off at bus stop 6. Figure 2B depicts the resulting clusters of drop-off points in the auxiliary graph, which has four clusters, one for the depot and one for each student.Each node number is provided inside the circle or square, whereas the number in parenthesis next to each node 1-9 indicates the node in Figure 2A that this node is a copy of.So, for example, nodes 5 in C 2 and 9 in C 3 in the auxiliary graph are copies of node 5 in the original graph, where students 2 and 3 can be dropped off.However, node 6 in the original graph is represented only by node 6 in the auxiliary graph, where only student 2 can be dropped off.

Model 3: Two-index formulation on the auxiliary graph based on MTZ
The following model is the two-index formulation on the auxiliary graph and the subtour eliminations are based on capacity MTZ inequalities.

Decision variables
x ij is a binary variable, which is equal to 1 if arc (i, j) ∈ A H is traversed; otherwise, it is zero.
i is a non-negative value that shows the arrival time to node i ∈ N H if node i is visited; otherwise, it is zero.
i is a non-negative value that shows the number of students on a vehicle immediately before visiting node i ∈ N H .If node i is not visited, then the value is undefined. min s.t.: The objective function (32) is to minimize the sum of arrival times, which contains the time that all students spend on the vehicles up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.Constraint (33) ensures that the students are dropped off at only one of the stops in their cluster.Constraint (34) ensures that |K| vehicles start travelling from the school.Constraint (35) ensures that a vehicle leaves stop i ∈ N H after visiting stop i.Constraint (36) illustrates the flow of the vehicles at each node i. Constraint (37) imposes the capacity control based on MTZ inequalities.Constraint (38) shows the arrival time consistency and imposes the connectivity for each vehicle.Constraints (39)-(42) show the domains of the decision variables.Now, by formulating the model as a GVRP, we can use the cluster-based MTZ constraints in our formulation that are defined first by [16], and then the stronger and lifted version of the MTZ capacity constraints for cluster-based formulations for the GVRP are introduced by [28].Therefore, instead of Constraint (37), we can take advantage of using the cluster-based formulations in our model.To show the cluster-based capacity constraint based on MTZ inequalities, let us first define u i as the number of students on the vehicle which visits cluster i, where i ∈ S.
For simplicity in demonstration, we define the following variables: Therefore, the lifted capacity constraints for the cluster-based formulations are as follows: It is now possible to replace Constraints (37) and (42) with Constraints (45)-(47) in Model 3. We note that Constraint (38) also eliminates subtours, except for when T ij = 0.In this case, inequalities (37) would serve to prohibit any subtours forming.
It is also possible to derive an alternative version of Constraint (38) which uses a variable t i as the arrival time to cluster i ∈ S instead of arrival time to a single node.

Proposition 2. The constraint
is a valid inequality for cluster-based time connectivity in the CSBRP.
Proof.If xij = 0, it means that there is not any connection between two clusters, then Tij = 0. Therefore, we will have the following valid inequality for any two clusters that are not connected.
When xij = 1, it means that there is a route connection between clusters i and j.Any cluster i may be on the first, intermediate, or last visit of a vehicle tour.If cluster i is the first visited cluster after the school, then the arrival time to cluster i is T0i , which shows the travel time from the school to the visited node inside cluster i.If cluster j is visited right after cluster i, then the arrival time to cluster j is the arrival time to cluster i and the travel time between the visited nodes of two clusters, which is equal to Tij .Therefore, ∀i, j ∈ S, Constraint (48) is valid.▪ Proposition 3. The constraint is a valid lifted version of Constraint (48).
Proof.We seek the largest value of  ji such that is valid.If xji = 0, this constraint is obviously satisfied for any value of  ji .If xji = 1, then xij = 0, and we obtain the following constraint: which provides the desired result as t j ≤ t i − Tji .▪ Proposition 4. The constraint is valid for the CSBRP.
Proof.Assume that P k is the length (in time) of the shortest path from the school up to node k ∈ N H .For any i ∈ S, the lower bound of t i is T0i , when a node inside cluster i is the first visited stop right after the school.In the case, when cluster i is not the first stop after the school, then t i is greater than or equal to the shortest path (in time) from the school to any visited node before cluster i and a travel time between the visited node and a node inside cluster i.Therefore, Constraint (53) can be obtained.▪ In model 3, we can substitute Constraint (38) with Constraints (50) and (53).Besides, we can use the value of T as explained in Proposition 1.

3.2.2
Model 4: Two-index formulation on the auxiliary graph based on GG The following is the two-index formulation on the auxiliary graph and the subtour eliminations are based on GG inequalities.This formulation is inspired by [11].

Decision variables
x ij is a binary variable, which is equal to 1 if arc (i, j) ∈ A H is traversed; otherwise, it is zero.
z i is a binary variable, which is equal to 1 when node i ∈ N H is visited; otherwise, it is zero.
f ij is a non-negative value, which shows the number of onboard students on arc (i, j) ∈ A H .If arc (i, j) is not traversed, then the value is zero. min The objective function ( 54) is to minimize the sum of arrival times of all students that contains the time that the students spend on the vehicle up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.In this objective, the total arrival time of the students is obtained by calculating the number of students on each arc.Constraint (55) shows that the students are dropped off only at one of the stops in their cluster.Constraint (56) shows |K| vehicles start travelling from the school.Constraint (57) ensures that when a vehicle visits stop i ∈ N H , it will depart from it.Constraint (58) shows a flow of vehicles at each node i. Constraints (59)-(62) impose the subtour elimination based on Gavish and Graves.Constraints (63)-(65) show the domains of the decision variables.

3.2.3
Model 5: Three-index formulation on the auxiliary graph based on MTZ In the following, we are going to introduce the three-index formulations on the auxiliary graph with the lifted version of cluster-based MTZ inequalities.

Decision variables
x k ij is a binary variable, which is equal to 1 if arc (i, j) ∈ A H is traversed by vehicle k; otherwise, it is zero.
The objective function ( 66) is to minimize the sum of arrival times for all students that contains the time that the students spend on the vehicle up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.Constraint (67) shows that the students will drop off in only one of the stops in their own cluster.Constraint (68) shows that |K| vehicles will depart from the school.Constraint (69) illustrates that after vehicle k visits node i, it will depart from the stop.Constraint (70) illustrates the flow of each vehicle at each node i. Constraints (71)-( 73) ensure that the maximum number of onboard students on each vehicle cannot exceed the vehicle capacity.Constraints ( 74) and (75) show the arrival time consistency and impose the connectivity for each vehicle.Constraints ( 76)-( 79) show the domains of the decision variables.

3.2.4
Model 6: Three-index formulation on the auxiliary graph based on GG The following is the three-index formulation on the auxiliary graph with GG subtour elimination constraints.

Decision variables
x k ij is a binary variable, which is equal to 1, if arc (i, j) ∈ A H is traversed by vehicle k; otherwise, it is zero.
is a non-negative value which indicates the number of onboard students, if vehicle k travels arc (i, j); otherwise, the value is zero.
The objective function ( 80) is to minimize the sum of arrival times for all students that contains the arrival time that the students spend on the vehicle up to the drop-off stops and the total time of walking from the drop-off stops to their home addresses.Constraint (81) shows that the students are dropped off in only one of the stops in their cluster.Constraint (82) shows |K| vehicles will depart the school.Constraint (83) illustrates that after vehicle k visits node i, it will depart from it.Constraint (84) illustrates the flow of each vehicle at each node i. Constraints (85)-( 89) show that the maximum number of on board students on each vehicle cannot exceed the vehicle capacity.Constraints (90)-( 92) show the domains of the decision variables.
In summary, an overview of the six proposed models is depicted in Figure 3.

COMPUTATIONAL RESULTS
In this part, we are going to show the computational results of the six mathematical formulations for the CSBRP.In Section 4.1, the data preparation is discussed, and the final results are provided in Section 4.2.

Data preparation
The proposed model was tested on real-life data from Innlandet in Norway.The dataset contains the locations of different schools, the students' home addresses, and the bus stops.To calculate the travel time between any pair of nodes, we used the Haversine distance metric, which is the distance between two points on the surface of a sphere measured along the surface of the sphere.We assumed that school buses have a speed of 30 (km/h) and the walking speed for the students is 5 (km/h).All models are applied to certain small-sized instances, which are for the schools where the number of students is in the set {7, 10, 13, 15, 16, 22}.By selecting five schools from the dataset that need to provide a service for the same number of students, we can evaluate the average performance of our proposed model across various school locations with equivalent student populations.We made the assumption that there were two vehicles with a capacity of 20 students available for all instances.The instances are accessible at https://github.com/Far-naz/SchoolBusRouting.
Here, the computational experiments on the set of data are described.All of the computational results were performed on a Windows platform running an Intel(R) Core(TM) i7-1260P @ 2.10 GHz laptop PC with 16 GB RAM.All the instances were solved by default settings of CPLEX (version 12.9), and the maximum running time is 1 h.

Computational results
In this section, we analyse the performance of the six mathematical formulations.Table 1 shows aggregated results for five different sets of schools.In Table 1, |S| is the number of students, W Max is the maximum allowed walking distance, |N| is the average number of nodes in models which are defined by the original graph, and |N H | is the average number of nodes for models defined by the auxiliary graph.Note that the total number of nodes in each instance refers to the number of potential stops, which primarily refers to students' home addresses and bus stops.An increase in the allowed walking distance causes a rise in the number of nodes (both |N| and |N H |), especially |N H | is larger than |N| when the allowed walking distance is greater than zero, because of copying the shared stops.
The other columns of Table 1 describe the obtained results for each model.#opt represents the number of instances (out of five schools) that reach optimality within the given time, CPU represents the average computation time (in seconds) for instances that reach optimality within 1 h, and if the total computation time for all five instances is over 1 h, then that value is displayed by "-." obj shows the average optimal or upper bound values (in hours).GAP is the average percentage of optimality gaps only for those instances that do not reach optimality within 1 h.Note that the optimality gap reflects the difference between the best known bound and the objective value of the best solution.For each instance, we indicate by bold values which model performs best, according to the following rule: first, we prioritize the average optimality gap, and as a second priority for tie-breaking we use the average CPU time.
All of the reported results are based on the introduced lifted version of the time connectivity and capacity constraints.

TABLE 1
Summary of results.Based on the obtained results, the maximum number of students that the models can solve to optimality is less than 22 students, and we stop running the models for larger instances.In general, among all the instances, Models 2 to 4 perform better than the other models.
Model 1 shows the worse performance even for very small-sized instances like seven students.As a result of the existence of big M in the formulation, this model is weaker in computation.Consequently, we stop running Model 1 for instances with larger than 10 students.
According to Table 1, it appears that the computation time of the models is strongly dependent on both the number of students and the maximum walking distance, such that by increasing these two criteria, the average computation time for all the models also increases.Based on the average computation times, Model 3 performs well only on very small instances like 7 and 10 students, while Model 4's performance is superior when there are more than 10 students.Models 3 and 4 are both based on a 2-index formulation.
On average, the number of instances that reach optimality is larger with Model 4 than with the other models.
Table 2 shows a brief summary of the computations grouped by instance size as given by the number of students.Column |S| is the number of students, Num is the total number of instances with that number of students, #opt is the total number of instances where optimality is reached, and CPU is the average computation time in seconds those instances where optimality is reached.
Table 3 shows details of the instances with 13 students.The columns of the tables are I which is the school index, N and N H are the number of nodes the and the auxiliary graph respectively, z * is the optimal or upper bound value, CPU shows the computational time ( in seconds) if the solution time is less than an hour GAP is the optimality gap ( in percentage).
To test the computation time of the models, we started running the models with W Max = 0, which means that all the students must be dropped off at their home addresses.Then we increased the value of W Max to check the obtained results and compare the performance of the models in obtaining the results with their computation times and the optimality gaps.As depicted in Table 3, by increasing the value of W Max , the computation time of the models is increased as well.We stopped running the models when the models cannot solve the problem within 1 h.
Table 3 shows how increasing the number of nodes cause increased computation times and optimality gaps.Especially the computation times of Models 5 and 6 are sensitive to the number of nodes.In comparing Model 3 and Model 5, which both are based on MTZ inequalities on the auxiliary graph, it is evident that the optimality gaps in Model 5 are larger than in Model 3.For instance, in school number 5 with W Max = 0.3, the reported gap for Model 3 is 16.39% and for Model 5 is 20.70%.In comparing Model 4 and Model 6, both of which are based on GG inequalities on the auxiliary graph, it is depicted that the optimality gap and computation time in Model 4 is smaller than in Model 6.When comparing Model 4 and Model 6 with Model 2, all of which are based on GG inequalities but differ in decision variables and graph definitions, Model 4 is superior to Model 2 both in terms of running time and optimality gaps.

Sensitivity analysis
For a sensitivity analysis, we run Model 4, as it performed better in comparison to other models.All of the instances in this section are based on the instances with 15 students in five schools.Figures 4 and 5 report the effect of the number of vehicles and the capacity of the vehicles on the computational time and the obtained objectives.
Figure 4A shows the effect of the number of vehicles (|K|) and the allowed walking distance (W Max ) on the computation time of Model 4. When only one vehicle is available, the average computation time is more than 1 h, which we did not report in the figure.However, it is evident from Figure 4A that by increasing the number of vehicles, the computation time is decreased.As a result, the model can provide a solution faster when there are more vehicles than when there is only one.Figure 4A also depicts that the computation time is also dependent on the value of W Max .The reason is that as the number of potential stops increases, selecting stops, and generating routes becomes more complicated for the model.

TABLE 3
The instance with 13 students.According to Figure 4B, computation times tend to increase as W Max increases.The figure also illustrates that when the vehicle capacity is limited, the computation time is increased.For instance, when W Max = 0.4 with two available vehicles, the computation time of the model to provide a solution with a vehicle's capacity of 8 passengers is 2500 (s) while with a capacity of 10 passengers is 1500 (s).
Figure 5 illustrates how the objective value varies in relation to the number of vehicles and their capacity and different walking distances.Based on Figure 5A, a larger number of vehicles available for transporting students reduces the sum of student arrival times.In Figure 5B, which is also based on two available vehicles, a smaller capacity of the vehicles increases the total arrival time of students.Due to the vehicles' limited capacity, the vehicles must travel longer route distances to deliver all the students to their drop-off points, resulting in longer total arrival times.

4.3.1
Analysing the CSBRP against the conventional SBRP In this subsection, we compare our proposed model with the conventional school bus routing problem.While the objective of our proposed model is to minimize the sum of arrival times of the students at their homes, the objective of the conventional SBRP is to minimize the total time or distance travelled by the buses.In order to compare the results, we calculated the sum of arrival times of the students from a solution obtained from the conventional SBRP problem, and call it converted total arrival times.Additionally, we calculated the total travel times by the school buses using a solution obtained from the CSBRP, and call it converted total travel times.
Figure 6 compares the average CPU time of the conventional SBRP and the CSBRP for 15 students in 5 different schools.It is depicted that the cumulative SBRP takes longer to run than the conventional SBRP, even though both models are based on two-index formulations with GG inequalities as subtour elimination, and only their objective functions are different.Figure 7A shows a comparison between the obtained results based on the CSBRP, and the converted total arrival times of the students from the SBRP.Based on Figure 7A, in general, the total arrival times in the SBRP are higher than the CSBRP, and the total arrival times increase with increasing allowed walking distances.
As shown in Figure 7B, the total walking time of the students is influenced by the value of W Max .In addition, the total walking time is greater in the SBRP in comparison to the CSBRP, as in the SBRP, the total travel time by the buses is more important than student convenience.
Figure 7A,B show that, in this example, the total arrival and walking times, respectively, for students remain constant in the CSBRP across varying values of W Max .Indeed, the route which visits all home addresses for W Max = 0 results also for the larger values of W Max , as it is faster for the students to be dropped off at the home addresses instead of having to walk home from another drop off location.However, as shown in Figure 7B, the total walking times in the SBRP increase as the allowable walking distance increases, revealing that in the conventional SBRP, the total length of the school bus is prioritized over the total arrival times of students to get to their homes.
Figures 8 and 9 show the obtained results from the CSBRP and the SBRP respectively.The results include the drop-off points of ten students with a maximum allowed walking distance of 0.4 (km) and the generated routes for two school buses with a capacity of ten passengers.
In Figure 8, the black square is the school and blue circles are the students' home addresses.This figure shows that all of the students are dropped off at their home addresses.There are two generated routes (in dashed blue and solid orange lines) where the total walking time of the students is zero, and the sum of arrival times for the students is 7.42 (h).
In Figure 9, the school is displayed by a black square, the students' home addresses are blue circles and the drop-off bus stops are red triangles.Besides, two routes from the model are displayed in Figure 9 by dashed blue and solid orange lines, where in the dashed line route, only student 8 is dropped off, and the other students are travelling on the solid line route.The number of bus stops for these ten students with W Max = 0.4 is 25, but in the figure only the visited drop-off stops are displayed.In this  solution, not all students are dropped off at their home addresses.Specifically, students number 2, 4, and 10 are dropped off at their home addresses, while student number 1 is dropped off at student 2's home address and the other students are dropped off at bus stops.The students who are dropped off at places other than their home addresses should walk to their homes.In Figure 9, the converted total arrival times for the students is 25.65 (h), the total walking time of the students is 0.35 (h), and the total travelled time is 6.53 (h).
In the SBRP, the two converted travel times are very different (6.422 and 0.113 [h]), while the two converted travel times obtained from the CSBRP are more similar (3.55 and 3.86 [h]) and as such avoid a long travel time for any individual student.
As explained earlier, the CSBRP involves a set of closed routes that start and end at the school.However, the objective only involves the total arrival times of the students, thereby excluding the traversal time of the edge that connects the last drop off students with the same home address would be serviced by different buses in an optimal solution to a unit-demand formulation, this solution would not be identified with general demands provided that demand splitting is not allowed.Moreover, if demand splitting were actually allowed in models with general demands, we consider that our unit-demand models are much more straightforward to implement.Second, it may well be that the real-world problem characteristics allow different walking times for different students with the same home address.This may generally be the case where rules related to the public transportation service specify different allowed walking times for students at different ages.So in this way, two or more students at different ages and with the same home address may well have different allowed walking times, which we can model by our unit-demand models but not by models with general demands.

CONCLUSION
In this article, we introduced a new objective for the school bus routing problem which has a focus on equity and fairness in the arrival time of students at their home addresses, when being transported by school buses from the school.We introduced six different mathematical formulations for the problem that differ in network design and decision variables.We also introduced lifted time constraints in the cluster-based formulations.Finally, the models' performance is compared on the set of real-life instances.We also compared our proposed model with the conventional school bus routing problem which has a focus on minimizing the total time of travel.The results show that in the conventional SBRP, students have to walk larger distances, and on average the students arrive later at their homes, while the objective of the CSBRP is to minimize the sum of arrival times of all students at their homes.
We can clearly observe that the mixed ILP approach in this article is only able to solve the smallest instances to optimality within the given time limit.Therefore, as further research, it can be investigated how to design exact and heuristic algorithms to solve the CSBRP for large-size instances.

FIGURE 2
FIGURE 2 Schematic view of the potential stops for three students on the original and auxiliary graphs in the CSBRP.In (B), numbers in parentheses refer to nodes in the original graph.(A) Subsets of nodes in the original graph.(B) Clusters of nodes in the auxiliary graph.

FIGURE 3
FIGURE3 Overview of the models.

FIGURE 5
FIGURE 5 Comparison between obtained objective based on number of vehicles and capacity of vehicles over different allowed walking distances.(A) Number of vehicles and allowed walking.(B) Capacities and allowed walking.

FIGURE 6
FIGURE6 CPU time differences between the conventional SBRP and the CSBRP.

FIGURE 7
FIGURE 7 Comparison between the obtained results from the SBRP and the CSBRP.(A) Total arrival times.(B) Total walking times.

FIGURE 8
FIGURE 8The obtained result from the CSBRP model for ten students.

FIGURE 9
FIGURE 9The obtained result from the SBRP model for ten students.

TABLE 2
Summary of computations grouped by instance size.
Comparison between CPU time based on number of vehicles and capacity of vehicles over different allowed walking distances.(A) Number of vehicles and allowed walking.(B) Capacities and allowed walking.